Universal scaling functions of critical Casimir forces obtained by Monte Carlo 
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Effective Casimir forces induced by thermal fluctuations in the vicinity of bulk critical points are 
studied by means of Monte Carlo simulations in three-dimensional systems for film geometries and 
within the experimentally relevant Ising and XY universality classes. Several surface universality 
classes of the confining surfaces are considered, some of which are relevant for recent experiments. 
A novel approach introduced previously [EPL 80, 60009 (2007)], based inter alia on an integration 
scheme of free energy differences, is utilized to compute the universal scaling functions of the critical 
Casimir forces in the critical range of temperatures above and below the bulk critical temperature. 
The resulting predictions are compared with corresponding experimental data for wetting fihns of 
fluids and with available theoretical results. 
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I. INTRODUCTION 

The confinement of a fluctuating medium generates ef- 
fective forces acting on the corresponding surfaces. Close 
to the critical point of a continuous phase transition the 
relevant fluctuating degree of freedom is the order pa- 
rameter of the phase transition. The effective force re- 
sulting from the confinement of such critical fiuctuations 
is known as the critical Casimir force fc ■ This force has 
a universal character in the sense that it is largely in- 
dependent of the microscopic details of the systems and 
of the confining surfaces but depends only on some of 
their gross features (which characterize the correspond- 
ing bulk and surface universality classes), as it is typically 
the case for bulk and surface critical phenomena. Such 
forces were first discussed by Fisher and de Gennes [l| 
on the basis of finite-size scaling 0] for a fluid system 
confined by two parallel walls. 

After early qualitative observations [1, 0] the first 
quantitative experimental evidence for such a force was 
provided by the study of wetting layers of "^He [Ej, where 
fc originates from the confined critical fiuctuations as- 
sociated with the superfiuid transition in the fluid film; 
fc adds to the omnipresent background dispersion forces 
which together determine the equilibrium thickness L of 
the wetting layers @ . The dependence of L on tempera- 
ture T provides an indirect measurement of fc', varying 
the undersaturation allows one to tune L and thus to 
probe the scaling properties of fc as function of T and 
L @, 0] • Later on, wetting layers of classical [1, Q and 
quantum binary liquid mixtures [loj have been studied 
and in two cases it has been possible to determine quan- 
titatively the critical Casimir force near a critical [8| and 
a tricritical [l0| point. Only recently, however, the exis- 
tence of the critical Casimir effect has been demonstrated 
by a direct measurement of the femto-Newton force be- 
tween a planar wall and a colloidal particle immersed in 



a near-critical binary liquid mixture 

The universality of the Casimir force fc allows one 
to investigate its temperature dependence via represen- 
tative models. Recently wc have briefly reported 
a novel approach for the Monte Carlo computation of 
the critical Casimir force which allowed us to study the 
scaling behavior of fc in the experimentally relevant 
cases mentioned above and to provide results for fea- 
tures of fc which were theoretically not accessible be- 
fore. Specifically, as follows from finite-size scaling the- 
ory [ij, the temperature and the geometry depen- 
dence of the critical Casimir force fc per unit area A 
and in units of ksT = (3^^ can be expressed in terms 
of a universal scaling function the form of which de- 
pends on the shape of the geometrical confinement, on 
the bulk universality class of the confined medium, and 
on the surface universality classes of the confining sur- 
faces The latter are related to the boundary condi- 
tions (BC) imposed by the surfaces on the rel- 
evant fiuctuating field, i.e., on the order parameter (OP) 
of the underlying second-order phase transition. 

Binary liquid mixtures near their demixing points 
belong to the bulk universality class of the three- 
dimensional (3D) Ising model, whereas liquid ''He near 
the superfiuid temperature of the critical end point of 
the A-line belongs to the bulk universality class of the 
XY model. In the aforementioned experiments involv- 
ing thin films of classical fiuids, both confining surfaces 
preferentially adsorb one or the other of the two com- 
ponents of the binary mixture. This corresponds to the 
surface universality class of symmetry-breaking surface 
fields [l^; the sign of the surface field (-1- or — ) acting 
at the boundary of the system indicates which compo- 
nent of the mixture is preferentially adsorbed. Accord- 
ingly, (H — ) BC reflect the fact that effectively the two 
surfaces attract different components of the liquid mix- 
ture, whereas (-1 — h) (and, equivalently, ( )) BC corre- 



spond to the case in which the two surfaces effectively 
attract the same component. In the case of the colloidal 
suspension studied in Ref. [HI both surfaces could be 
treated chemically such that (++) as well as (H — ) BC 
have been realized. For the wetting experiment of Ref. Q 
the appropriate BC are (H — ). In the case of wetting ex- 
periments for pure superfluid ^He the superfluid OP 
vanishes at both interfaces; there are no surface fields 
which couple to the superfluid OP. This corresponds to 
the symmetric Dirichlet-Dirichlet BC {O, O) based on the 
so-called ordinary (O) surface universality class. 

Due to the complexity of technical challenges as well 
as due to conceptual issues like the dimensional crossover 
in three-dimensional films, theoretical studies of the scal- 
ing functions of the critical Casimir forces by analytic 
means have been either limited to mean-field calcula- 
tions or have been confined to the disordered phase or to 
BC without symmetry breaking fields. Therefore Monte 
Carlo simulations offer a highly welcome tool to overcome 
these shortcomings and to study, inter alia, the aforemen- 
tioned experimentally relevant universality classes within 
the whole temperature range. 

Our computer simulations of the critical Casimir force 
are based on the integration scheme of free energy differ- 
ences via the so-called "coupling parameter approach" 
and we computed the scaling functions for the 3D Ising 
model with (-l-f), (-1 — ), Dirichlet-Dirichlet (O, O), and 
periodic BC (PBC), as well as for the 3D XY model with 
Dirichlet-Dirichlet (O, O) and periodic BC. In all cases 
we studied the film geometry. The experimental data of 
Refs. and @ turn out to be in a good agreement with 
our simulation results which are, in addition, consistent 
with those obtained by alternative numerical approaches 
based on the computation of either the expectation value 
of a suitable lattice stress tensor [l^ for the 3D Ising 
model with periodic BC, or of the internal energy den- 
sity, followed by an integration over the temperature [l^l , 
for the XY model with (O, O) BC. We also find good 
agreement with the results of the de Gennes-Fisher local- 
functional method extended to the Ising universality with 
(++) BC 

The purpose of the present study is to elucidate the 
relevant details of the approach used in Ref. and to 
present new results for both the Ising and the XY bulk 
universality class in three dimensions. In particular, we 
extensively discuss the important issue of corrections to 
scaling and the fitting procedure necessary to obtain the 
estimates of the scaling functions d from the raw MC 
data. Several functional forms of corrections to scaling 
are considered and the ensuing differences in the result- 
ing scaling functions are described. In particular, the 
estimates for the universal Casimir amplitudes at Tc are 
obtained. 

Our presentation is organized as follows: In Sec. HIl we 
provide the basic theoretical background, i.e., the mod- 
els, the critical Casimir force, and the scaling functions 
are defined. In Sec. IIIII we summarize our method for 
the computation of the scaling functions. New data for 



the XY model with (O, O) as well as with periodic BC 
are presented in Subsec. IIV Al They have been obtained 
for larger lattices and with a better accuracy compared 
to the results presented in Ref. Discussions of the 

dependence of the corresponding scaling functions on the 
aspect ratio of the simulation cell and of the corrections 
to scaling are included. For the case of (O, O) BC in the 
XY model we present the comparison with the experi- 
mental data for wetting films of ^He [5 [and with the MC 
simulation results obtained in Refs. 3 Data for 

periodic BC in the XY model are compared to the avail- 
able field-theoretical predictions above the bulk critical 
temperature Tr [1, 0, [13, HH and to the MC simulation 
data of Ref. [1^1 . The analysis of the 3D Ising model is 
reported in Subsec. IIV Bl where we present new data for 
the Casimir scaling function for the (O, O) BC, the as- 
pect ratio dependence of the Casimir scaling functions for 
periodic BC, the determination of the universal Casimir 
amplitude via the analysis of the finite-size corrections, 
and the detailed description of the fitting procedure. In 
addition we compare our results for periodic BC in the 
Ising model with recent field-theoretical predictions for 
the behavior of the corresponding scaling function above 
Tc @, 0, HOj HH and with results in two dimensions (2D) 
For (-1— 1-) and (-1 — ) BC we provide a comparison of our 
data with the exact results in 2D [22| and with mean-field 
predictions [2^ as well as with results of the extended 
de Gennes-Fisher local-functional method applied to the 
case of (-t--l-) BC [l^, [l^. The experimental data for the 
scaling function obtained from the wetting experiments 
for a binary liquid mixture in Ref. Q are compared with 
our MC results for (-1 — ) BC. We end with a summary 
and conclusions in Sec. IVl 



II. THEORETICAL BACKGROUND 

We consider the Ising and the XY model defined on a 
three-dimensional simple cubic lattice via the Hamilto- 



3j , 



(1) 



where J > is the spin-spin coupling constant, the sum 
(i, j) runs over all nearest neighbor pairs of sites i and j 
on the lattice. In the Ising model, has only one com- 
ponent Si e { + 1,-1}, whereas in the XY model is 
a two-component vector with modulus |si| = 1. Tem- 
peratures and energies are measured in units of J. The 
inverse critical temperature is (3c ~ 0.2216544(3) [2^ for 
the Ising model, whereas pc = 0.45420(2) Q for the 
XY model. We consider film geometries, i.e., lattice cells 
of sizes Lx X Ly x Lz with Lx = Ly ^ Lz = L and 
A ~ Lx X Ly, with periodic BC in the x and y directions 
(in which the system has linear extensions and Ly). 
In the z direction we consider (O, O) and periodic BC for 
the XY model and fixed, (O, O), and periodic BC for the 
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Ising model. The (++) and (H — ) BC are realized by fix- 
ing the boundary spins to values Si = +1 (+) or = — 1 
(— ) whereas (O, O) BC arc realized by free surface spins. 

In a film geometry with thickness L and large trans- 
verse area A, the Casimir force /c per unit area A and 
in units of hgT = is defined as 



fciP,L)^-dr/dL, 



(2) 



where ^""{(3,1) = (3L[f - /''""'(/3)] is the excess free 
energy which depends on the type of the BC, / is the 
free energy of the film per volume V = LA and /^"""^ 
is the bulk free energy density. From the general theory 
of finite-size scaling Q and based on renormalization- 
group analyses 0, Q we expect the Casimir force to take 
the universal scaling form 



(3) 



where the scaling function depends on the spatial 
dimension d and on the BC. Here r = (/3c — f3)/P = 
(T — Tc)/Tc is the reduced temperature and f = ^^ItI""^ 
is the bulk correlation length which controls the spatial 
exponential decay of the two-point correlation function. 
The critical exponent i/ equals 0.6301(4) and 0.662(7) 
for the Ising and the XY bulk universality class in three 
dimensions, respectively [1^ ; are nonuniversal ampli- 
tudes above {+) and below (~) with = 0.501(2) gi| 
for the Ising model on the simple cubic lattice, whereas 
^(J" = 0.498(2) m for the XY model. The values of ^+ 
quoted here refer to the amplitude of the second moment 
correlation length ^2'"^l however, £,/S,2^d — 1 for (3 < (3c 
for both the Ising and the XY model [2i,[2i|. 

At T = Tc the scaling function reduces to the uni- 
versal Casimir amplitude ^ (0) = {d — 1)A, which has 
been extensively studied in the literature (see, e.g., 
Refs. @, S [H, [H, [H, iOl). Determining the whole 
temperature dependence of the scaling function and its 
dependence on the spatial dimension d is a much more 
challenging task. 

For the Ising universality class with (O, O), (++), and 
(H — ) BC in the film geometry theoretical results are 
available in d = 2 from the exact diagonalization of 
the transfer matrix [22| and in d > 4 from mean-field 
theory [2^. In d = 3 theoretical results are available 
for T > Tc and periodic BC investigated both by MC 
simulations (at Tc) [l^ and by field-theoretical meth- 
ods SSSEII as well as for Dirichlet [1,0], von Neu- 
mann BC [0, 0l7and Robin BC [13] investigated by field- 
theoretical methods. Recently, the extended de Gennes- 
Fisher local-functional method has been applied in order 
to study the case of (++) BC within the full temperature 
range 

For the bulk universality class of the XY model in 
film geometry with (O, O) BC theoretical results for the 
Casimir force scaling function are available in d = 3. 
They include field-theoretical calculations for tempera- 
tures T > Tc \E , [3 and numerical results from MC simu- 
lations (l^ . [l7l |T In the low temperature limit the specific 



features of the superfluid ^Hc were taken into account 
in Ref. [2^ and the contribution to the Casimir force 
resulting from the capillary-wave like fluctuations on the 
surface of ''He wetting films was determined. For T < Tc, 
which corresponds to temperatures below the superfiuid- 
normal fiuid transition temperature Tx of the A tran- 
sition, certain qualitative features of the Casimir scaling 
function have been recently understood within the frame- 
work of the Landau- Ginzburg mean-field theory [29l.[30j. 

For large areas A, the total free energy F{(3, L, A) of 
the confined system can be written as 

F{f3,L,A) EE ALf = A[Lf^^'\/3)+r'r{P,L)]. (4) 

The quantity contains two L-independent surface 
contributions in addition to the finite-size contribution 
f°^{P,L) — /'^'^(/3,oo) the L-dependence of which gives 
rise to the effective Casimir force. On a lattice ('), the 
derivative in Eq. ^ is replaced by a finite difference and 
fc{(3,L) is given by 



+pr'\P), (5) 



where AF(/3, L, A) = F{f3, L, A) - F{l3, L - I, A). One 
can consider different definitions of the lattice derivative 
than the one we have implemented in Eq. Different 
choices give rise to different corrections to the leading 
behavior of the Casimir force scaling function. 



III. METHOD 
A. Computation of free energy differences 

Monte Carlo methods are generally not efficient for 
the computation of quantities, such as the free energy F, 
which cannot be expressed as ensemble averages. Nev- 
ertheless, free energy differences, such as AF{(3, L, A) 
we are interested in, can be cast in such a form via 
the so-called "coupling parameter approach" (see, e.g., 
Ref. [3l|). This is a viable alternative to the method 
used in Ref. [ll] in which a suitable lattice stress tensor 
has been defined in such a way that its ensemble average 
renders AF. So far, however, this latter method is only 
applicable for periodic BC. 

If one is interested in the Monte Carlo computation 
of the difference Fi — Fq between the free energies Fi ~ 
— ^ In exp(— /3/^i) (i £ {0,1}) of two lattice models 
Mo and Mi with the same configuration space C but 
different Hamiltonians Hq and Hi , respectively, it is con- 
venient to introduce an "interpolating" system A4cr(A) 
with the crossover Hamiltonian 



i/cr(A) = (l-A)Ho + Ai/l, 



(6) 



where A G [0,1], and again the same configuration 
space C. As a function of the coupling parameter 
A, i?cr(A) interpolates between Hq and Hi as A in- 
creases from to 1 and accordingly the free energy 
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a) 



FIG. 1: Bond arrangement for the computation of the free 
energy difference in Eq. Q (see main text). The crossover 
Hamiltonian Her (c) belongs to a system which interpolates 
between those described by the Hamiltonians Hq (a) and Hi 
(b). 



i^cr(A) = -ilnX;cCxp(-/37?er(A)) of X„(A) interpo- 
lates between and Fi. The difference Fi — Fq can be 
trivially expressed as Fi — Fq = dAi^^'^(A) where 
is the derivative of Fcr(A) with respect to the coupling 
parameter: 



dX 



Xcr(A) 



(7) 



which takes the form of the canonical ensemble average 
(. . .)a<„(a) of AH = Hi — Hq and therefore it can be 
efficiently computed via MC simulations of the lattice 
model A^ci (A). As a result one can conveniently express 
the difference in free energies as an integral over canonical 
averages (see, e.g., Ref. [3l[): 



Fi — Fq 



A\ (AH) 



Mor(A) 



(8) 



According to Eq. the Casimir force wc arc inter- 
ested in is related to the difference AF{/3, L, A) between 
the free energies F{(3,L,A) and F{P,L — 1, A) of the 
same model on two lattices with different numbers of 
sites and therefore different configuration spaces. In or- 
der to apply the method described above for the compu- 
tation of Ai^(/3, L,A) one identifies the model M.q^ its 
Hamiltonian, and the associated configuration space C 
with the corresponding ones of the model we are inter- 
ested in on the lattice A x L, as depicted in Fig. [IJa), 
so that Fo(/3,L, A) = F(/3,L,A). The final system Mi 
has to be chosen such that it has the same configura- 
tion space C as M.q. This is achieved by adding to 
the model on the lattice A v. (L — 1) - for which we 
want to compute the free energy F(/3, L—\,A) ~ a two- 
dimensional lattice of size A with suitable degrees of free- 
dom and lateral periodic BC (see Fig. [IJb)). The Hamil- 
tonian Hi of A^i is defined such that the added layer 
does not interact with the remaining part of the system 



and therefore Fi(/3, L, A) = F{f3, L-1,A) + F2d(/3, A), 
where F2d(/3,A) is the free energy of the isolated two- 
dimensional layer. This layer can be thought of as the 
one at position ko G {1, 2, . . . ,L} (along the z-direction) 
in the model Mq which then decouples from the rest of 
the lattice upon passing from A = to A = 1, i.e., from 
Fig. [1] (a) to (b). The resulting crossover Hamiltonian 
Hcr{X) (see Eq. ([5])) additionally depends on the origi- 
nal position ko of the extracted layer. In particular, in 
the three-dimensional models we are mainly interested 
in, the fluctuating degrees of freedom are one- (Ising) or 
two-component (XY) vectors s^^y^z — where i = (x, y, z) 
specifies the lattice site — which interact only with their 
nearest neighbors on the same lattice, with a coupling 
strength J ~ \ (indicated by solid bonds in Figs. [T] (a) 
and (b); J is absorbed into (3). For them one explicitly 
finds 



Ai? = Hi — Hq — — 2^(Sx,y,fco-l ■ Sx.y.fco + l 

^x^y^ki:) — ! ' ^x,y,kQ ^x,y,kQ ' ^x,y.kQ + l^ • 

(9) 

The resulting iJcr(A) = Hq + XAH is characterized by 
the coupling constants depicted in Fig. [TJc). The free 
energy difference AF (see Eqs. ([5]) and ([5])) can be finally 
expressed as 

AF(/3, L,A)^~I{P,L,A) + F2U A) (10) 

where I{l3,L,A) = dX{AH) ^^^(^x)- Note that 
H„{X) (see Fig. [BJc)), AH (see Eq. ®), and therefore 
{AH)h„{x) depend on the value of fco whereas I{/3, L, A) 
is actually independent of it, as long as the boundary 
conditions are not affected by the extraction of the ko-th 
layer as A varies between and 1. For fixed and open 
BC in the z-direction this requires ko 1, whereas for 
PBC there is no such a restriction on ko and (AH) fj^^i^x) 
is actually independent of it. In our simulations we have 
chosen fco ~ L/2. 

Once AF(/?, L, A) has been computed, one has still to 
subtract /*'""'(/?) from it (see Eq. in order to de- 

termine the Casimir force in a film of assigned thickness 
L — 1/2. However, the accurate computation of the bulk 
free energy density /'^""^(/S) is a numerical problem by it- 
self and extracting it from finite-size data requires a very 
accurate analysis. In order to avoid this complication in 
the computation of the Casimir force, it is convenient to 
consider the difference between the forces acting in slabs 
of thicknesses Li and L2 > Li. 

Afcil3,Li,L2,A) 

^/c(/3,^i-^,A)-/c(/3,L2-i,A) (11) 



[I{P,Li,A)~Iif3,L2,A)] 



in which the contributions of both /''""^(/S) and -F2d(/3, ^) 
actually cancel. Accordingly, the procedure to calcu- 
late the scaling function of the Casimir force consists 
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of the following steps: (1) For a given geometry L x A 
and temperature via MC simulations we compute 
the ensemble averages {AH)h^,.{x) for different values of 
A G {Xi,...Xn}- (2) On the basis of these N values 
we calculate the integral I{(3, L, A) in Eq. ([8]) via numer- 
ical integration. (3) These computations are repeated 
for different sizes L, A, and temperatures f3~^, yielding 
numerical estimates for A/c(/3, Li, L2, A) (see Eq. (jlip ). 
(4) The scaling function -d in Eq. ^ is retrieved from 
the numerical data for A/c as described below. The 
results presented in Sec. IIVI have been obtained by us- 
ing the Simpson integration method with = 20 men- 
tioned above in step (2) and by using pairs of geometries 
(Li, L2) = (i, 2L) with L = 13, 16, 20 for the Ising model 
and L = 10, 15, 20 for the XY model, as introduced in 
step (4) and for fixed aspect ratios p = L/^/A. (The 
motivation for our choice L2 = 2L will be provided in 
the following subsection.) The method of Ref. [l3| takes 
advantage of the possibly available numerical knowledge 
of the bulk energy density u^^^^ of the model of interest 
whereas here the analogous information on the bulk free 
energy f^^^^ is not required for the determination of the 
Casimir scaling function, making our approach applica- 
ble also to cases in which there is no detailed knowledge 
of u*"""^ and /'^""^. 



B. Determination of the scaling function 

The scaling function -d of the Casimir force can 
be extracted from the temperature dependence of 
A/c(/3, Li, L2, A), for fixed Li^2 and A, by using the fact 
that fc in Eq. (fTT|) scales according to Eq. ^ for large 
Li^2 and A. In order to highlight these scaling properties 
it is convenient to introduce the quantity 

g{y;L,,L2,A)^{L,-l/2fx ^^^^ 
AfcW^ P{y;L,),Li,L2,A), 

as a function of y, where l3{y;Li) = (3c/[l + y{Li — 
112)-^/"]. According to Eq. © and with r = (/3c-/3)//3, 
g is expected to scale as 

g{y- L,,L2, A) = e{y) - a-''9{a'/''y) , (13) 

where a — {L2 — l/2)/(Li — 1/2) is the width ratio, and 
§{y) is the Monte Carlo estimate of 0{y) = 'd{y/{£,'^Y'''); 
here d = 3. Note that, even though is independent of 
this geometrical realization of the simulation cell, 9 might 
depend on it via A and Li_2 due to corrections to scaling. 
For a given pair of geometries LixA and L2XA, the avail- 
able Monte Carlo data for AfciP, Li, L2, A) at differ- 
ent temperatures allow one to determine g{y; Li, L2, A) 
for a discrete set of values of y. In order to deter- 
mine 9(y) from the numerical data for g{y; Li, L2, A) 
with fixed L1.2 and A, one can solve Eq. (fT3|) itera- 
tively. One can expect (see below) that this yields a 



solution for L2 > Li, i.e., a > 1 together with the prop- 
erty 9{\y\ 00) (which holds apart from T < Tc 
in the XY model). As a first approximation of the ac- 
tual 9{y) one takes 9o{y) = g{y; Li, L2, A), which can 
be improved by taking into account that Eq. (jl3p yields 
9{y) ^ 9o{y) + a-''9{a^/''y) ^ 9o{y) + a-''9o{a^/''y). Ac- 
cordingiy, a better approximant 9i{y) is provided by 

9i{y)^9o{y)+a-''9o{a^^''y). (14) 

The values of at the point for which no MC 

data might be available, are obtained by cubic spline in- 
terpolation of the available ones. In Eq. one can 
replace 9o by using Eq. ITSll . yielding 9i{y) — 9{y) — 
a-^'^9{a^/''y) ~ 9{y) - a-^'^9i{a'^/''y), which indicates 
how the approximant 9i (y) can be improved by introduc- 
i^gUv) = 0i(y)+a-2rf^i(a2/^'y) = 0(y)-a-4d0(a4/.y), 
This expression can in turn be used to further improve 
the approximant along the same lines. The resulting it- 
erative procedure yields a sequence of approximants 

4>i(y) = 0k-i{y) + a-^''"%-i{a^''"^''y), (15) 

which converges very rapidly because the correction to 
the fc-tli approximant is of the order of i.e., ex- 

ponentially small in 2^^ and, in addition, 9{y) is gener- 
ally expected to decay exponentially for large |y|. With 
a ~ 2, already for fc = 5 one has a^^ ~ 3.5 x 10~^^ in 
three dimensions (d = 3). The choice of a ~ ^2/^1 is a 
compromise between two competing aims: a small value 
reduces the sizes of the lattices required for the computa- 
tion of /c but on the other hand it decreases the accuracy 
of a given approximant in determining 9. With our choice 
of geometries (Li,L2) — {L,2L), one has a ~ 2 and a 
very good approximation of 9{y) = 9k^oo{y) is already 
provided by 9z{y). 

C. Details of the MC simulations and test of the 
method 

In order to compute the canonical average {AH)^^^^x) 
we use a hybrid MC method which is a suitable mixture 
of Wolff and Metropolis algorithms . Specifically, for 
the Ising model each hybrid MC step consists of four flips 
of a Wolff cluster according to the Wolff algorithm, typ- 
ically followed by 3A attempts to flip a spin Sx,y,z with 
z € {ko — 1, ko, ko + 1 1, which are accepted according 
to the Metropolis rate [34|- An analogous method, with 
a suitable implementation of Metropolis and Wolff algo- 
rithms, has been used for the XY model [l^l, i.e., a flip of 
a Wolff cluster according to the Wolff algorithm is typi- 
cally followed by the implementation of moves according 
to the Metropolis algorithm [3^. 

In order to test the program we have computed numer- 
ically g{y, 5, 10, 9) as a function of y for the Ising model 
on a lattice 3 x 3 x L with periodic, (-f-f ), and (-) — ) BC, 
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finding perfect agreement with the result of the analytic 
calculation based on the transfer-matrix method. 



D. Corrections to scaling 

Finite-size scaling is known to be valid asymptoti- 
cally for large lattices and small values of r, i.e., a 
large correlation length ^ Q • Away from the asymptotic 
regime corrections to the leading (universal) scaling be- 
havior become relevant. These non- universal corrections 
affect both the scaling variables and the scaling func- 
tions and depend on the details of the model as well as 
on the geometry and the boundary conditions [H, [33 |. 
Renormalization-group analyses reveal that there is a 
whole variety of sources of corrections which arise from 
bulk, surface, and finite-size effects 0]. 

For the limited thicknesses L of the lattices we inves- 
tigated with our MC simulations it is necessary to take 
corrections to scaling into account in order to obtain data 
collapse [13, [13 ■ III the present case, the finite-size scal- 
ing variable t[L / ^'^Y^^ (in the following associated with 
the reduced temperature t) is expected to acquire a lead- 
ing correction of the form 

x^T{L/^+f/^{l + g^L--), (16) 

where uj is the leading correction-to-scaling exponent in 
the bulk which takes the values 0.84(4) and 0.79(2) ^ 
for the three-dimensional Ising and XY universality class, 
respectively. Corrections to the scaling behavior of the 
critical Casimir force /c are expected to be of the form 

/c(/3,L,A)=L-^^(x,L-') 

~L-^t9(.T)[l + L-">(x) + ...], 

for L 1, where the exponent lo' controls the leading 
corrections to the scaling behavior of the lattice estimate 
/c. Its value is determined by that irrelevant surface or 
bulk perturbation of the Hamiltonian H which has the 
smallest scaling dimension and which also affects /c. In 
the generic bulk case one has uj' = oj. But its value can be 
suitably increased (so that the infiuence of the corrections 
is reduced) by using improved Hamiltonians and observ- 
ables, which can also serve as representatives of the same 
universality class. This is described in detail in Ref. [1^. 
In the presence of surfaces, irrelevant surface perturba- 
tions might yield w' < w, but we are not aware of either 
theoretical or numerical studies of this issue. In addition, 
for small lattice sizes, next-to-leading corrections to scal- 
ing might also be of relevance. If uj' > 1/2, these correc- 
tions are generically provided by analytic terms ^ 
(even though they might be absent in some quantities). 
The interplay between the leading and next-to-leading 
corrections (especially if they are sizable) might result in 
an effective exponent ujcs- The current accuracy of our 
Monte Carlo data and the relatively small range of sizes 
L investigated here do not allow a reliable determination 



of Lo' and (j^ix). In particular it will turn out that the cor- 
rections to scaling are quite well captured by assuming 
(/){x) ~ f/2, i-e., a constant within the range of the scaling 
variable we have investigated, and an effective exponent 
Woff for the size dependence. 

In the discussion of the expected scaling behavior of 
/c we have assumed that the aspect ratio p = Lj \J~A is 
small enough (i.e., p <C 1) so that the scaling behavior 
in Eq. ([3]) holds, which formally corresponds to the limit 
A CO. On the other hand, the actual Monte Carlo sim- 
ulations have been performed on lattices with small but 
non-zero p and therefore possible additional, p-dependent 
corrections have to be taken into account in order to 
be able to extrapolate our results to the limit p ^ 0. 
The numerical results in Ref. [sst on the (universal) p- 
dependence of the Casimir amplitude i?(0, p) (see Eq. ([3])) 
of the three-dimensional XY model with periodic and free 
boundary conditions suggest ?9(0, 0) ~ '&{Q, p)(l+rp'^) for 
p < 0.5, where r is a constant. (This is confirmed also by 
the analysis in Ref. [l3| . ) In what follows we assume that 
this dependence on p carries over to the whole scaling 
function so that i?(a;, 0) ~ '0{x, p){l + (j)2{x)p'^). Although 
the amplitude 02 (a^) of the correction might depend on 
the scaling variable x (and possibly on L^'^ ), we shall as- 
sume that 4'2{x) ~ r2, i.e., a constant at least within the 
range of values of the scaling variable x which is studied 
in the present analysis. On the same footing, we expect 
a quadratic p-dependence of the finite-size scaling vari- 
able x{p) ~ a;(0)(l 4- rip^) associated with the reduced 
temperature t, where ri is a constant and a;(0) is given 
by Eq. (fTB|) . Taking into account all these corrections, 
we identify 

x = T(^Ay (i + g^L-")(l+rip2), (18) 

as the finite-size scaling variable, in terms of which the 
expected scaling behavior of /c is given by 

/c(/3, L, A) = L-\l + g2L-"-0(l + r2p')-H{x). (19) 

We shall aim at fixing the non-universal constants ri.2, 
g^j, and 52, which generally depend on the boundary con- 
ditions, in such a way that the data collapse of the avail- 
able Monte Carlo data is optimal. 

In most of the cases considered below, the accuracy of 
the data and the range of sizes L investigated do not allow 
for the reliable determination of both the amplitude (72 
and the exponent ujcs of the correction. Therefore we fix 
tJoff = \ — lo' , which actually leads to a reasonably good 
data collapse within the considered range of the scaling 
variable. 

In the absence of corresponding dedicated theoretical 
and numerical analyses, there is no a priori reason why 
one should prefer the use of a specific form of correc- 
tions to scaling, because all of them amount to an ef- 
fective way of accounting for these corrections. Accord- 
ingly, adopting a pragmatic approach, we shall choose 
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that form which leads to the best data collapse or to 
the best fit. Specifically, we use the following functional 
forms of corrections to scaling: 

case (i): 

/c(/3,L, A) - L-\l + giL-')-\l + r^p^Y^S (x) , 

(20) 



case (ii): 



case (iii): 



(21) 



(22) 



and 



case (iv): 
fc{l3,L,A) 



Case (i), (ii), and (iv) become all equivalent for large 
lattice sizes L. On the other hand, for smaller lattice 
sizes, they lead do different estimates. The coefficients 
5ii52,<73 and 31,32 are determined in such a way as to 
optimize the data collapse in the resulting estimate for 
'd{x) (see below). The factor of the form (iv), with two 
fitting parameters, will be considered only if data cor- 
responding to several different values of L are available, 
so that the resulting estimates for gi and 32 are reliable. 
Case (iii) of corrections to scaling works well for the XY 
and the Ising model with periodic BC. In cases in which 
corrections to scaling are not small, the ansatz used for 
their dependence on L might lead to a biased estimate of 
the scaling function •d{x). 

In order to highlight and assess the relevance of the 
different kinds of corrections, we present in the following 
sections also the MC data for the function g{y; L, 2L, A) 
which is the primary quantity determined by our MC 
simulation and from which the scaling function 'd(x) is 
eventually obtained according to the procedure described 
in Subscc. fill Bl In the absence of corrections to scaling, 
data for g (sec Eq. dl^ ) with L2 = 2Li; Li = L, as a 
fimction of y = {(3^/13 - 1)(L - 1/2)^/" = t{L - 1/2)^^" 
with fixed p = L/^/A but different sizes L should col- 
lapse on a single master curve, which, however, it is not 
always the case (see, c.f.. Fig. HJa)). In order to account 
for the corrections to scaling we proceed as follows: First, 
for fixed values of L and A = {L/ pY we determine the 
Monte Carlo data for g (see Eq. for different val- 

ues of the inverse temperature (3. Second, from the plot 
of g as a function of the rescaled reduced temperature 
y = t(L — 1/2)^/", i.e., from g{y; L, 2L, A), we determine 
the estimate of the scaling function 9{y), according to the 
procedure described in Subsec. IIII Bl This procedure is 
repeated for the different geometries considered in each 
case. Because of corrections to scaling and corrections 



due to p 7^ 0, the resulting estimates 0{y) actually de- 
pend on the specific values of L and A = {L/p)'^, i.e., 
6 = 6{y] L, p). In order to extract the asymptotic limit 6 
of the scaling function of the Casimir force from the lat- 
tice estimate 6, we account for corrections in accordance 
with Eqs. (HH]) and (with possibly different forms for 
the L-dependent corrections, see Eqs. which 
involve several fitting parameters. In those cases in which 
we apply corrections to scaling due to the aspect ratio de- 
pendence of the function gijj] L, 2L, A) (which turns out 
to be the case only for the XY model), the actual fitting 
procedure we shall use is divided into two steps. 



In the first step we fix the value of L (L = 10 for the 
XY model) and consider data corresponding to different 
aspect ratios p {p^^ = 4, 5, 6, 8, 10 for XY). The parame- 
ters ri and r2 are therefore determined such that the data 
for (1 + r2p^)~^9{y; L, p) (cx 9{x) for fixed L) as a func- 
tion of 2/(1 -I- rip^) (cx X, see Eq. (fT6|) . for fixed L) yield 
the best data collapse onto a single curve which ideally 
corresponds to the scaling function in the limit A 00, 
but which is still affected by L-dependent corrections to 
scaling. (This procedure actually assumes that, accord- 
ing to Eqs. (fT5|) and ri and r2 do not depend on 
L.) 



In the second step we fix the value of p (p = 1/6 for 
both XY and Ising) and we determine g^ and 32 (or gi 
or both gi and 32- depending on the specific form as- 
sumed for the corrections) in such a way that the data 
for (1 + g2L^'^''")0{y; L, p) (oc 6{x) for fixed p) as a func- 
tion of y(l -|- g^L~'^) (cx x, see Eq. (fTB)) . for fixed p) 
yield the best data collapse onto a single curve. (This 
procedure actually assumes that, according to Eqs. (fT8|) 
and (Uni), in a first approximation corrections to scaling 
do not depend on p.) The details of the fitting proce- 
dure are described in Appendix [X] Our final numeri- 
cal estimate of the scaling function 'd{x) of the Casimir 
force is then provided by the curve which result from 
plotting (1 -I- s(2i"'^'''0(l + 'r2P^)~^Q[y]L,p) (or equiva- 
lent forms as given by the cases (i)-(iv)) as a function of 
y{l + g^L-'^){l + rip2)/(^+)i/'^ = x, where the fitting 
parameters have been fixed according to the procedure 
described above. 



Finally it is worthwhile to keep in mind that besides 
the common corrections to the leading critical behavior 
in experimental data, the available experimental results 
for critical Casimir forces contain an additional source 
of corrections in that the thickness L of the (wetting) 
films is definition-dependent up to a microscopic length 
if) [3^. Accordingly, only the leading term is univer- 
sal whereas the correction ~ (-o/L is even definition- 
dependent. Moreover, also the relation between the ex- 
perimental values iexp and the theoretical values Lthco 
suffers from the same kind of uncertainty. 
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IV. RESULTS 

In this section wc summarize the numerical results for 
the scaling function of the critical Casimir force within 
the three-dimensional XY fSubsec. lIV Ap and Ising (Sub- 
sec. IIVBP universality classes with different boundary 
conditions. As mentioned in the Introduction, the former 
are relevant for the interpretation of the experiments with 
wetting films of ^He [H , whereas the latter apply to the 
case of classical binary mixtures [1, . In most of the 
presented plots the size of the symbols are of the order 
of the statistical error. In these cases the corresponding 
error bars are not shown in the figures. 



A. XY model 

For the simulations of the XY model we have consid- 
ered films of thicknesses L = 10, 15, and 20, and trans- 
verse areas A = x Ly ~ 6L x 6L corresponding to an 
aspect ratio p = L/\/A ^ 1/6. At the boundaries in the 
X- and y-directions we impose periodic BC, whereas in 
the z-dircction we consider cither free surface spins, cor- 
responding to the (O, O) universality class, or periodic 
BC. 

In Fig. [2] we report the data corresponding to 
g{y;L,2L,A) (see Eq. ^) for (a) (0,0) and (b) pe- 
riodic BC. Corrections to scaling, which are signaled by 
the fact that data corresponding to different L do not fall 
onto the same master curve, are much more pronounced 
for the case of (O, O) BC (see Fig. EJa)) as compared 
with the case of periodic BC (see Fig. [IJb)). The same 
holds for the dependence of the data on the aspect ratio 
p (data for {O, O) arc not shown, data for periodic BC 
are presented in Fig.[3|). For (O, O) and periodic BC cor- 
rections to scaling arc more relevant for y < j/min, where 
j/min is the value of y at which the function g{y] L, 2L, A) 



attains its minimum. For y 



> 



Vn 



the data obtained 



for different L follow a common curve. We note that the 
bulk correlation length ^ of the XY model is infinite for 
all temperatures below Tc- £,{T < Tc) = oo. Accord- 
ingly, within the XY model the scaling variable y can be 
expressed as y = [{L - l/2)£_+ f^]^^" only for T > T^. 

Interestingly, for both types of BC the aspect ratio de- 
pendence is particularly strong in the range of tempera- 
tures around the minimum of the function gly; L, 2L, A), 
i.e., —2 < y < —1 and — 1 < < for (0,0) and peri- 
odic BC, respectively (see Fig. [3]). According to Fig. [2l 
the minimum of g{y; L, 2L, A) for periodic BC occurs at 
the reduced temperature Tmin = —yimn/{L — ^Y^'^ — 
-0.31/(L- ^y^". For (0,0) BC it occurs slightly fur- 
ther away from Tc, i.e., at ymin — —1-34 and —1.50(3), 
depending on the value of L. We find that changing L 
at a fixed aspect ratio results in slight relative shifts of 
the data whereas changing p at fixed L leads to much 
more pronounced differences (see Fig. This behavior 
is expected to be related to the finite-size effects near 





^ -0.2 
0.4 



(N 





0^ 



-0.6 
-0.8 
-1 



L=lO \ 
L=15 \ 
L=20\ ' 
XY ^" ,\ 
(0,0) BC 



(a) 



-6 -5 -4 -3 -2-10 1 

y 





-0.1 
-0.2 
-0.3 
■0.4 



^ -0-5 
-0.6 



XY 

periodic BC 



1 / L=10 \ 
(b) V L--20 1 



-4-2 2 

y 



FIG. 2: Monte Carlo data for g{y = r(L - i)^/"; L, 2L, ^ = 
(L/pf) (see Eq. JH}, r = {T - Tc)/Tc) in the three- 
dimensional XY model for L — 10, 15, 20, and fixed inverse 
aspect ratio p^^ = 6. In (a) and (b) we present the result 
for (0,0) and periodic BC, respectively. For T > Tc, i.e., 
y = y+ > one has y+ = [{L - ^)^o/^+V^''- ^^r the XY 
model ^_ = cxD for all temperatures T < Tc- The Kosterlitz- 
Thouless transition of the two-dimensional film occurs at 
y = yc,oo = -2.69(3) and y = yc,P = -0.996(1) H 
in (a) and (b), respectively. 



the thin film critical point. Within the Ising model, 
for an infinitely large transverse area A the point at 
which the film with {O, O) or periodic BC exhibits the 
2D critical behavior is located on the bulk coexistence 
line _ff = at a size-dependent temperature Tc(L) < 
such that ^(T = Tc{L)) ~ L. Accordingly, upon in- 
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FIG. 3; Monte Carlo data for g{y = r(L - i)^/"; L, 2L, yl = 
(L/pf) (see Eq. g!]), t = {T ~ Tc)/Tc) within the three- 
dimensional XY model with periodic BC for L = 10 and dif- 
ferent values of the inverse aspect ratio p^^ . For y > one has 
y = y+ = [(L-i)Cj/^+]i/''. For the XY model 5- = oo for all 
temperatures T < Tc. The Kosterlitz-Thouless transition of 
the two-dimensional film occurs at j/ = yc,p = —0.996(1) [4l| . 
Note the enlarged scales as compared with Fig. [2jb) . 



creasing L, Tc{L) approaches the three-dimensional bulk 
value Tc as Tc(L ^ oo) = T^il + ycL-^^") 0], where 
Uc is negative, non-universal, and depends, inter alia, 
on the type of BC. The corresponding scaling variable 
x{L) = [{TciL) - Tc)/Tc]{L/i+f/'' tends to a universal 
and BC-dependent value x* = x{L co) = Vciia )^'^^'' ■ 
Hence, yc is expected to lie in the vicinity of the mini- 
mum of the function (?(?/; L, 2L, A). Accordingly, around 
its minimum the function g{y; L, 2L, A) should exhibit a 
strong dependence on the aspect ratio p if the bulk cor- 
relation length S^^^ , associated with the shifted critical 
point of the two-dimensional film [2^ , becomes compa- 
rable with the characteristic transverse length = a/A 
of the simulated system. 

Within the XY model the critical point of the thin 
film belongs to the Kosterlitz-Thouless (KT) universality 
class d^j. The KT theory predicts that upon approach- 
ing this critical point from the high temperature phase 
the correlation length £^kt ^ exp[(l — (3/f3^'^)~'^ ], 
i'^^-^ = 1/2, diverges exponentially. The shift of Tc{L) 
relative to the bulk critical point is expected to scale 
with the film thickness L in the same way as for the 
Ising model, i.e., {Tc{L) - Tj/Tc ^ ycL^'^/" for large 
L, with 1/ = 0.662(7) for the 3D XY model. This pre- 
diction is in agreement with MC simulations of vari- 
ous models belonging to the XY universality class and 
confined in films with free [H, [H, H^l or periodic [4lj 
BC. The critical exponent v obtained from early simu- 
lations 3^. 39| of films with free BC was slightly larger 
(i/ ~ 0.7 isTsgj) than the theoretically predicted value 



ly = 0.662(7), due to rather strong corrections to scal- 
ing which need to be taken into account in order to ob- 
serve the theoretically expected behavior The re- 
sults of Ref. [13 for Xqq = —7.64(15) yield the estimate 
y^ OO = -(Cc}")i/''a;5o = -2.69(3) for the location of the 
shifted KT transition in the XY model with free bound- 
ary conditions, whereas yc,p = —0.996(1) for PBC I4ll . 
It turns out that in the simulations with free BC [38[ 
the positions of the maxima of the thermodynamic func- 
tions such as the peak of the specific heat or the peak 
of the susceptibility do not coincide with the transition 
point but occur at ca. 1.3 x Tc{L). (These quantities are 
not related to singularities of XY films.) With increas- 
ing film thickness L the absolute distance in temperature 
of these peaks from Tc{L) decreases and for L = 10 the 
simulations of Ref. [s^ report a shift of less than 10%. 
There is also experimental evidence that as a function of 
temperature the position of the minimum of the Casimir 
force of '*Hc films, which belong to the universality class 
of XY films, coincides with the position T,„(L) of their 
specific heat maximum (see Subsec. VD and Figs. 21 
and 32 in Ref. whereas the onset of superfluidity 
in these films occurs at Tc{L) < T^iL) (see Subsec. VD 
and Figs. 24, 32, and 33 in Ref. [4^). A similar behavior 
may be expected to hold for the function g{y; L,2L, A). 
Indeed, as can be seen from Figs. [5] and [31 the minima of 
the function g{y; L,2L, A) lie in the vicinity of the cor- 
responding values of yc- Therefore, similar to the Ising 
model, the strong aspect ratio dependence around the 
minimum might occur when the exponentially diverging 
bulk correlation length , associated with the KT criti- 
cal point of the film, becomes comparable with the char- 
acteristic transverse length Ly = ^/A of the simulated 
system. 



1. Dirichlet-Dirichlet boundary conditions 

We consider first the case of (O, O) BC. As evidenced 
by Fig. [2lja), in order to achieve a good data collapse 
of the curves corresponding to different lattice sizes we 
have to account for corrections to scaling according to 
Eqs. (flS)) and (fT9)) . As a phenomenological ansatz for 
the effective corrections we take Woff = 1 and consider two 
functional forms for the i-dependent corrections to the 
scaling function: case (i) [Eq. ([^ ] and case (ii) [Eq. (PTjl ] 
as discussed in Subsec. IIIIDI As a result of the fitting 
procedure, in the interval x € [-6,-2.1] (see Eq. (US])) 
we find ri = 1.18(10), = 2.40(13), gi = 5.83(25), and 
g^ = 2.25(15) in case (i) and g2 — —2.98(8) in case (ii) 
with the same values for ri , r2 , and guj as in case (i) . Fig- 
ure [4] shows the corresponding resulting estimates of the 
scaling function dix) of the critical Casimir force. The 
quality of the data collapse for the two cases separately 
clearly indicates that Eqs. ((20)) . ((2T|) . and HH) are very 
effective ways of accounting for the corrections to scaling 
in this system. We find that {^{x) is slightly affected by 
the choice of the functional form of corrections to scaling 
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and indeed in the two cases one finds estimates of d{x) 
which have the same shape but the overall amplitude is 
reduced by a factor R ~ 0.9 in case (ii) as compared with 
case (i). The dashed line represents the scaling function 
which has been determined in Ref. [l3| on the basis of a 
different numerical method and assuming corrections to 
scaling of the form (i). Even though this result is actually 
biased by that particular choice (a point which has not 
been discussed in Ref. [l3l ) , the very good agreement be- 
tween the different approaches provides a highly welcome 
independent test of both methods. 

Our MC results for d{x) compare well also with the 
experimental data of Ref. [5[ . (For a meaningful compar- 
ison between the numerical and the experimental scal- 
ing function, the abscissa ri^/" of the experimental 
data presented in Ref. Q has to be properly normal- 
ized as T(L/^(j''°'^^^)^/'' by using the experimental value 
^+(exp) ^ 1 432A [43, [43.) In particular, the position 
of the pronounced minimum of the scaling function is 
properly captured. The corrections to scaling of form (i) 

yield xl^l = -5.43(2) and ^ 79(x«J = -1.396(6), 
whereas those of form (ii) result in x[^^^^ ~ —5.43(2) 
and ^tZ = ^(^min) = -1-260(5). The corresponding 
experimental values are a;|^j^jj|''' = —5.7(5) and = 
— 1.30(3). Taking into account the aforementioned bias 
affecting the results of Ref. [l3l and the sensitivity of 
the resulting scaling function to the assumed form of the 
corrections to scaling we conclude that our estimates for 
Xniin and iJmin are compatible also with those presented 
there (—5.3(1) and —1.35(3), respectively). As expected, 
due to the presence of the Goldstone modes below Tc, 
both the experimental and the MC data do not approach 
zero for x — > —00 but saturate at some finite negative 
value at low temperatures. However, the absolute value 
of the saturation as obtained from the MC simulations 
is smaller than the experimental one. This difference, 
which extends deep into the non-critical regime, is, inter 
alia, due to ^He specific properties and to the occurrence 
of capillary waves on the liquid- vapor interface of the crit- 
ical ^He wetting films. This point has been discussed in 
Ref. [11]. In Fig.|4]the gray vertical bar indicates the uni- 
versal value Xq Q — —7.64(15) of the scaling variable cor- 
responding to the occurrence of the Kostcrlitz-Thoulcss 
transition at T = Tc{L) in the film, as inferred from 
MC simulations of lattice models in the XY universal- 
ity class presented in Ref. The Kostcrlitz-Thoulcss 
transition is accompanied by an actually invisible essen- 
tial singularity ~ exp(— const/ yj\x — x*\) in the behavior 
of the specific heat which, as discussed above, displays 
a pronounced maximum at a temperature T = Tm{L). 
Accordingly, one does not expect to find any particular 
signature of this transition in the scaling function of the 
Casimir force for a; ~ a;*, in distinction to the case of the 
Ising model (c.f., Subsecs. l!VB2l and UVBSl) . 

Finally, for completeness, in Fig. |4] we have also in- 
cluded (dash-dotted line) our mean field result for the 
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FIG. 4; Scaling function {too of the Casimir force for the 
three-dimensional XY model with (O, O) BC. The MC data 
reported in this figure refer to lattices with L = 10, 15, and 
20, with fixed inverse aspect ratio 1/p = 6. Corrections to 
scaling have been accounted for according to two different 
ansatze, provided by Eq. (|20fl and Eq. (|21|l : the correspond- 
ing numerical results are denoted by (i) and (ii) , respectively. 
With corrections to scaling of the form (ii) , the shape of the 
resulting scaling function is almost indistinguishable from the 
one obtained with corrections to scaling of the form (i), but 
its overall amplitude is reduced by a factor ii ~ 0.9. For (i) 
our MC data compare very well with the corresponding ex- 
perimental data from Ref. (solid line) and with the MC 
data of Ref. [l3] (dashed). Due to the Goldstone modes 
1^00(2; —00) — const 7^ 0. The dash-dotted line shows the 
mean field scaling function [29l . [30l | normalized to the depth of 
the minimum of the MC data (i). The levelling off of the ex- 
perimental data for x — » —00 contains a component which 
is specific for *He wetting films [2^ and cannot be captured 
by an XY lattice model. The gray bar indicates the posi- 
tion and uncertainty of the universal value Xq q = —7.64(15) 
of the scaling variable x corresponding to the occurrence of 
the Kosterlitz-Thouless transition in the film, as inferred from 
MC simulations of lattice models in the XY universality class 
presented in Ref. [40||. 



Casimir scaling function '^^q^"^^ obtained from the lim- 
iting case of the vectoralized Blume-Emery-Griffiths lat- 
tice model corresponding to the model of pure ^He [29j . 
The scaling function is normalized to the depth of the 
minimum of the MC data. For large L, ^^00"^^ agrees 
very well with the ones obtained from the 0{2) Landau- 
Ginzburg continuum theory [2^, [s^l . 



2. Periodic boundary conditions 

In this subsection we discuss the XY model with pe- 
riodic BC. According to Fig. [IJb) corrections to scaling 
are much less pronounced in this case than for (O, O) 
BC (Fig. [2Ia)), suggesting that the exponent ujcs might 
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FIG. 5: Critical Casimir amplitude Ap{L) for the three- 
dimensional XY bulk universality class and periodic BC, esti- 
mated from lattices of several thicknesses L and inverse aspect 
ratio 1/p = 6. Due to corrections to scaling, Ap depends 
on L. The solid line represent the best fit to the numer- 
ical data based on Eq. (|24|) and allows one to extrapolate 
the value of Ap(L) to the scaling limit L ^ oo, resulting 
in Ap = —0.2993(7) (■). With • we indicate the numerical 
estimate Ap = —0.28 provided in Ref. [l^ . 



be actually larger than 1. In addition, the dependence 
of the numerical data on the aspect ratio p turns out 
to be relevant only in the restricted range —^SuS^ 
of the scaling variable (see Fig. [3]), so that the assumed 
forms of the aspect ratio corrections in Eqs. and 
do not work best. In the present case, the accuracy of 
our Monte Carlo data allows us to study in some detail 
also the Casimir amplitude A = i?(0)/2. Upon focus- 
ing on such a quantity in a broader range of geometries 
(6 < L < 20) it turns out that for this amplitude the 
corrections to scaling are not properly accounted for by 
the previous ansatze (case (i) and case (ii), Eqs. ([20|) and 
(PT|) . respectively). We have therefore tried also a fit of 
the exponent ujcs according to Eq. p9|) with ri^2 = 
(case (iii), Eq. (|22|)). which yields for the Casimir ampli- 
tude 



A{L) = A{1 + gsL- 



(24) 



With this ansatz, our data for Ap{L) are very well fit- 
ted for uJcs = 2.59(4) and = 14.9(7) in the interval 
< L~-^ < 0.15. (At present, the origin of this rather 
large value of tJcff is not clear.) The comparison between 
the numerical data and the fit is reported in Fig. [5l The 
value of the Casimir amplitude extrapolated to the scal- 
ing limit L — > oo is Ap(oo) = Ap = —0.2993(7) which is 
slightly smaller than the previous estimate Ap = —0.28 
(see Ref. [lB| and the discussion below). Note, however, 
that our estimate is biased by the particular form Eq. ((24|) 
assumed for the corrections to scaling. 



The analysis of the Casimir amplitude Ap{L) suggests 
that the corrections to scaling for periodic BC are well 
captured (in the range of sizes and of the scaling variable 
investigated here) by Eq. ([^ (case (iii)) and Eq. 
with ri = 0. The resulting estimate for the scaling func- 
tion 'dp is reported in Fig.[6]for which we adopt the values 
for ^3 and locS which we determined from the analysis of 
the correction to scaling for Ap{L). It turns out that 
a very good data collapse is achieved even without cor- 
recting the abscissa, i.e., with ~ 0, ri ~ 0, within 
the range of the scaling variable x we have investigated, 
which actually includes the interval — 1 < y < in which 
the corresponding function g shows a more pronounced 
dependence on p. 

As another valuable test of the method, our results are 
compared with the corresponding MC simulation data 
obtained previously in Ref. within a different ap- 
proach, i.e., by computing the average value of the lat- 
tice stress-tensor. In Fig. |6]we report the data set corre- 
sponding to the lattice size L = 20 investigated therein. 
The shapes of the two scaling functions are very similar 
but the data points from Ref. [lB| are shifted upwards 
with respect to the ones we have obtained. This discrep- 
ancy might be due to the uncertainty in the normaliza- 
tion factor used in Ref. [ll|, where the vertical scale of 
the data for dp has to be adjusted on the basis of an 
independent estimate. This estimate has been obtained 
from the e = 4 — d-expansion of the ratio Ap„/Ap i 
of the Casimir amplitudes for 0{n) models with the re- 
sult Ap,„=2 = -0.28 so that ^?p(0) = 2Ap,„=2 = -0.56. 
In contrast; the method presented here provides abso- 
lute values of the amplitude and the scaling function. 
In addition to the uncertainty concerning the normaliza- 
tion factor, in Ref. [1^1 no corrections to scaling have 
been applied in the determination of i?p. The present 
MC results provide the estimates Xmin = —0.73(1) and 
'&p{xynin) = —0.633(1) characterizing the position of the 
minimum of the scaling function. 

For the scaling function of the XY model with peri- 
odic BC some analytical predictions are also available; 
for a thorough comparison of the scaling function ob- 
tained within various approaches see Ref. [4^ . Here we 
discuss only the comparison with the recent results based 
on a suitable perturbation theory for the 0{n) model in a 
film geometry withperiodic BC [13, [2l| , which improves 
previous analyses p, Q of this scaling function for T >Tc 
by taking into account a higher-order contribution to the 
perturbation theory which involves fractional powers of e. 
In the case n = 2 (XY model) and in agreement with our 
MC data this latter analytically available scaling func- 
tion decreases monotonically for a; — > and thus allows 
for the formation of a minimum below Tc (without being 
able to reach it) whereas the previously available ana- 
lytic scaling function exhibits a minimum above Tc- The 
analytically estimated value for the critical Casimir am- 
plitude is Ap ~ -0.43 (i.e., dp{Q) ~ -0.86) which in 
absolute value is larger than the MC result. In Fig. [S] 
this analytically predicted scaling function is reported, 
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for comparison, as a solid line. 

As already mentioned above, one characteristic feature 
of the scaling function of the critical Casimir force in the 
XY model (and, more generally, in systems with con- 
tinuous symmetry) is its saturation at a nonzero nega- 
tive value d^x — > oo) < at low temperatures, which 
occurs for all non-symmetry breaking BC. This is due 
to the fact that, even well below the critical temper- 
ature, the fluctuations of the order parameter exhibit 
long-ranged correlations due to the Goldstone modes as- 
sociated with the broken continuous symmetry, which re- 
sult in a non-vanishing long-ranged Casimir force. For 
periodic BC the saturation value i?p(— oo) is signifi- 
cantly more negative and is approached more rapidly 
than in the case of {0,0) BC. The line of arguments 
presented in Ref. [1^ for the theoretical calculation (TH) 
of i?[^o^(-oo) -C(3)/(87r) ~ -0.049 (disregarding ad- 
ditional helium-specific surface fluctuations) can be ex- 
tended to the present case by considering periodic (in- 
stead of Neumann as in Ref. [2g) BC for the fluctuations 
of the phase field of the order parameter in the film. In 
three dimensions this yields 



i?^'"\-oo) = 2A^''^ ~ -0.38, (25) 



where Ap = — C(3)/(27r) ~ —0.19 is the Casimir ampli- 
tude for a one-component {N — 1) fluctuating Gaussian 
field in a film with PBC (see, e.g., Eq. (9.2) in Ref. 01), 
so that dp ' {—oo)/'dQ q' {—oo) = 8. The numerical data 
corresponding to the MC simulations presented in Fig. [6] 
yield '&p{—oo) = —0.383(4) (obtained by fitting the data 
points in the region — 14 < a; < —10, two of which are 
actually not shown in Fig. [6l with a constant). This 
is in very good agreement with the theoretical predic- 
tion ^"'(-oo) in Eq. Note that the MC data of 
Ref. [3 give -dp^—oo) ~ —0.33, a value which is biased 
by the choice of the normalization of the scaling function, 
as mentioned before. We point out, however, that the line 
of arguments in Ref. [1^ assumes that, deep in the low- 
temperature phase, the phase field obeys Neumann BC 
and that the magnitude of the complex order parameter 
(superfluid density) is spatially constant across the film, 
i.e., that the effects of the surfaces are effectively negli- 
gible. This might not be the case in the presence of the 
Goldstone modes which can cause the magnitude of the 
order parameter vary algebraically within the film. 

Finally, in Fig. [6] wc report as a gray vertical line 
the universal value x*p ~ —2.82(2) of the scaling vari- 
able X corresponding to the occurrence of the Kosterlitz- 
Thouless transition in the film, as inferred from the MC 
simulations of the XY model in a film with PBC [4l| . As 
in the case of (O, O) BC, there is no singularity possibly 
visible in dp{x) associated with this transition. 




FIG. 6; Scaling function i9p of the Casimir force for the three- 
dimensional XY model with periodic BC. The corrections to 
scaling are taken into account by Eq. (|22p (case (iii)) and 
Eq. (|18|) with ri = 0. The shape of our MC data compares 
very well with the corresponding MC data (•) of Ref. [l^ . For 
a discussion of the relative shift of the data sets see the main 
text. The solid line corresponds to the analytical prediction 
in Ref. [2l[. Due to the Goldstone modes, in agreement with 
Eq. ([251), J?p(a; — > -oo) = -0.383(4), see horizontal dashed 
line. Contrary to (O, O) BC in Fig. |1 for periodic BC MFT 
yields 'd^^^'^\x) = for a; ^ 0. The gray vertical line in- 
dicates the position of the universal value x*p = —2.82(2) of 
the scaling variable x corresponding to the occurrence of the 
Kosterlitz-Thouless transition in the film, as inferred from 
MC simulations [iH ]. 

B. Ising model 

In the case of the Ising model we have determined the 
scaling function i? for (H — ), (-I--I-), Dirichlet-Dirichlet 
(O, O), and periodic BC. The first two BC are relevant for 
interpreting the results of the experiments in Refs. [l,[ll| 
which use as a critical medium classical binary liquid mix- 
tures near their demixing point. 

In our simulations we have used lattices with L = 10, 
13, 16, and 20 and with = Ly = 6L, i.e., p = 1/6. 
Each data point has been averaged over at least 10'^ hy- 
brid MC steps. 

1. (++) and (H — ) boundary conditions 

We first discuss the cases of (4--t-) and (H — ) BC, for 
which we find that in the critical regime the numerical 
data for the function g{y; L,2L, A) are practically inde- 
pendent of the aspect ratio p = Ll\fA (see Fig. [7]). The 
presented data correspond to L = 10 and to inverse as- 
pect ratios p"^ =6, 10, 14. In the case of (H — ) BC the 
aspect ratio becomes relevant for y < —4, where the be- 
havior of the system is dominated by the presence of 
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FIG. 7: Plot of g{y;L,2L,A = (L/pf) (see Eq. ^) for the 
three-dimensional Ising model with L — 10 and 1/p = 6, 10, 
and 14. (a) and (b) refer to (+- 1-) and (H — ) BC, respectively, 
and the coincidence of data points corresponding to differ- 
ent values of p demonstrates that the geometry of the lattice 
does not affect the resulting finite-size critical behavior in the 
region —4 < y < 10. 



the strongly fluctuating interface which separates the re- 
gions with predominantly positive and negative magne- 
tization. The extent of these fluctuations is known to 
be particularly sensitive to the spatial extension and to 
the geometry of the system in the directions parallel to 
the interface (i.e., in the and Ly directions); there- 
fore the aspect ratio p plays an important role for these 
fluctuations. In d = 3 one expects a strongly increasing 
parallel correlation length ^|| which governs the decay 
of the correlations in the direction parallel to the inter- 
face, i.e., ^ exp(L||/(4^)) with Ly = = Ly 
In addition, these strong interfacial fluctuations cause 

the scaling function to decay to zero for x —oo 

much more slowly than the scaling function (see, 



c.f.. Figs. [To] and 11). 

Contrary to the aspect ratio, L-dependent corrections 
to scaling are rather important for the Ising model with 
(-I--I-) and (H — ) BC. By using the phenomenological 
ansatze in Eqs. (HH]) and (HOI) or ([HI) with ri,2 = (which 
account for the negligible dependence of the data on p) 
we have obtained a good data collapse for the scaling 
functions calculated for L = 13, 16, and 20. However, 
these ansatze fail to describe the data for the critical 
Casimir amplitude A in the broader range of thicknesses 
6 < L < 20, as it is the case of the XY model with pe- 
riodic BC. It turns out that the corrections to scaling in 
this range are very well captured by the functional de- 
pendence (iv) (Eq. introduced in Subsec. HlIDi with 
''1,2 = 0) which for the critical Casimir amplitude yields 



A(i) = A 



(1+ 32^-1) ■ 



(26) 



As in the case of the XY model with periodic BC we 
shall determine the parameters gi and g2 (according to 
Subsec. IIIIPp for both (+-I-) and (H — ) BC on the basis 
of the analysis of the corrections to scaling to the cor- 
responding Casimir amplitude and then these values are 
employed in order to calculate the scaling functions 
and t9_| 

In Figure [5] we present numerical data for A as a func- 
tion of 1/L for both {++) (a) and (+-) (b) BC with the 
corresponding fit carried out according to Eq. (|26|) in the 
interval < 1/L < 0.1. Other variants of the fit function 
(within the same fit interval), such as A(l -I- giL~^)~^ 
and A{1 + g2L~^), indicated as (i) and (ii), respectively, 
are also presented for comparison. 

For (++) BC the fitting parameters are gi = 
—2.6(1.2), 52 = 6.6(3.7) and the resulting estimate for 
the Casimir amplitude is A+-|_(L s- oo) = A-|_+ = 
-0.376(29), i.e., t?++(0) = -0.75(6), which compares 
quite well with the previous MC result '!?++(0) = 
-0.690(32) [H shown as a full circle in Fig. [H^a); for 
the latter result corrections to scaling were not taken 
into account. Field-theoretical predictions 1^4.+ (0) = 
—0.652 ... — 0.346 give numbers slightly smaller in ab- 
solute value which depend on the approximant used to 
re-sum the field-theoretical e — A ~ d-expansion up to 
0(e) series (see Ref. for detafls). 

For (-) — ) BC we have found gi = —1.8(1), .92 = 

8.54(43) and we estimate A_| (L 00) 

2.71(2), i.e. 



(0) = 5.42(4), in agreement with the 

y^j - 6(2) [I but slightly 

larger compared to the previous MC estimate {f^ (0) = 

4.900(64) [2II (indicated as a fuU circle in Fig.lSfb), stiU 
affected by finite-size corrections) and the analytical esti- 
mates 'd^+]:\0) = 3.16 . . .4.78. The latter depend on the 
approximant used to re-sum the 0{e) series (see Ref. [2^ 
for details). 

By using the values of gi and 52 obtained previously in 
the context of the Casimir amplitude we determine the 
coefficient g^ of the correction to the scaling variable x 



experimental value i}^^^\0) 
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FIG. 8: MC data for the critical Casimir amplitude A{L) 
of the three-dimensional Ising model with (a) (++) and (b) 
(H — ) BC, as a function of the inverse lattice size L (for lattices 
with fixed inverse aspect ratio 1/p — 6). L-dependent correc- 
tions to scaling give rise to the dependence A(L) such that 
A+_|./+_ = A_|.+/_|._(L — > oo). The solid line corresponds to 
the best fit obtained by using the fitting ansatz in Eq. p6|) in 
the interval < l/L < 0.1. For comparison we present also 
the best fits using the ansatze (i) A{L) = A{l + giL~^)~^ and 
(ii) A(L) = A(l -I- g2L~^). Our estimates (■) for the asymp- 
totic values of the Casimir amplitudes compare reasonably 
well with previous MC results (•) from Ref. [2j and with re- 
sults (♦) obtained from the de Gennes-Fisher local-functional 

(see Eq. (|T8|) with ri = 0) in order to achieve a good 
data collapse for the whole scaling function, with the re- 
sults = 2.04(15) for {++) BC and = 2.90(15) for 
(-1 — ) BC. The comparison between three phenomenolog- 
ical ansatze for the corrections to scaling, i.e., cases (i) 
[Eq. (ii) [Eq. dH])], and (iv) [Eq. are presented 

in Figs. M and [TUl for (++) and (-1 — ) BC, respectively. 
The scaling functions corresponding to the rational ex- 



pression for the corrections to scaling ansatz (case (iv)) 
lie in between the two others. 

Currently, for the film geometry with (++) BC there 
are no experimental data available for comparison, but 
in Fig. [5] can be compared with the prediction of 
mean-field theory [2^ (MET, solid line, normalized such 

that 47^)(0) = <"^.^'(0) [= 2Af+^^ see Eig. Ha)]) 
and with the prediction of the two-dimensional Ising 
model HI] (dashed line). Recently, the de Gennes-Eisher 
local-functional method has been extended to study the 
three-dimensional case with (++) BC [H, [H. In this 
latter (non-perturbative) approach one takes advantage 
of the knowledge of the values of bulk critical exponents 
and amplitude ratios in order to fix completely certain 
parameters of an effective model which is then used 
to calculate the structural properties and the free en- 
ergy of the system first in the presence of a single wall 
and eventually in thin films, giving access to the scal- 
ing function for {++) BC. The resulting scaling func- 
tion (dash-dotted line in Eig. [5]) is in very good agree- 
ment with the one (bottom set of data points in Eig. [5]) 
determined numerically via MC simulations by assum- 
ing corrections to scaling of the form given by Eqs. 
and HH) with ri.2 = and suitable values for the fit- 
ting parameters g^^ and gi (see above). This agreement 
suggests that corrections to scaling are properly cap- 
tured by such ansatze even for L — > oo. The predic- 
tion of the de Gennes-Fisher local-functional method for 
the critical Casimir amplitudes (shown as diamonds in 
Figs, li^a) and (b)) is A++ = -0.42(8) and A+_ = 3.1 
[l8| . which compares quite well with our MC results for 

A++ = —0.376(29); whereas for A^ the agreement with 

our result A^ = 2.71(2) is slightly less good. 

In the case of (H — ) BC we can compare with the 

experimental results of Ref. [Ij, with the prediction of 
mean-field theory [23| and with the corresponding result 
for the two-dimensional Ising model [l^ (see Fig. [TO]). 
The solid line, normalized similarly as for (-I--I-) BC, rep- 
resents the ME result, whereas the dashed line refers to 
the two-dimensional Ising model [1^ . We expect the ex- 
perimental data in Ref. @| to be affected by corrections 
to scaling already for x > 2, due to the relatively small 
corresponding value of ^/i < 30, where £ ~ 3A is the 
molecular scale set by the specific binary liquid mixture 
used in Ref. Q. In view of these difficulties, the compari- 
son between the MC and the experimental data in Eig. [TO] 
can be regarded to provide an encouraging agreement. 

Within the Derjaguin approximation our numerical re- 
sults for and form the basis for the calcula- 
tion (llj of the corresponding scaling functions for the 
critical Casimir potentials in the sphere-plate geometry, 
which turn out to be in remarkably good agreement with 
the actual experimental results for that geometrical set- 

ting [m. 

Comparing the scaling functions for d = 2, 3, and 4 
(MET) one finds that in the case of {++) [(+-)] BC the 
position of the minimum [maximum] moves away from 
the bulk critical point a; = as the spatial dimension in- 
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FIG. 9; Scaling function '&++ of the critical Casimir force in 
the three-dimensional Ising model with (++) BC and zero 
bulk field. Data points refer to lattices with fixed inverse 
aspect ratio 1/p = 6. The bottom and top data sets have been 
obtained by accounting for corrections to scaling according to 
Eq. (f20)) (case (i)) and Eq. ([21} (case (ii)), respectively. The 
intermediate data set, instead, considers corrections of the 
rational form given by Eq. (|23p (case (iv)). In each case the 
data collapse turns out to be very good within the range of the 
scaling variable x covered in the figure. The final estimate of 
the scaling function is biased by the functional form assumed 
for the corrections to scaling. The position Xmin — 5.90(8) of 
the minimum is insensitive with respect to these choices for 
the form of the corrections. For comparison we provide the 
prediction of mean- field theory [SJ (solid line), normalized 
such that i?^^^'^>(0) = [Figj8l;a)], the exact result 

for the two-dimensional Ising model [22| (dashed line), and the 
result from the extended de Gennes-Fisher local-functional 
method [l9l | (dash-dotted line). Note that the actual phase 
transition of the film occurs at a nonzero value of the bulk 
field. 



FIG. 10: Scaling function '&+- of the critical Casimir force 
in the three-dimensional Ising model with (H — ) BC and 
zero bulk field. Data points refer to lattices with fixed in- 
verse aspect ratio 1/p = 6. For comparison we provide the 
mean-field prediction [23| (solid line), normalized such that 
^(MFT)|,Q-| ^ ^(_MC)|,Q^ 2A+_, Fig.lSTb)], the exact result 
for the two-dimensional Ising model [22| (dashed line), and 
the set of experimental data points from Ref. 0|. The top 
and bottom data sets have been obtained by accounting for 
corrections to scaling according to Eq. (|20p (case (i)) and 
Eq. (|21|l (case (ii)), respectively. The intermediate data set, 
instead, considers corrections of the rational form given by 
Eq. (|23|l (case (iv)). In each case the data collapse turns out 
to be very good for x > —20. The final estimate of the scal- 
ing function is biased by the functional form assumed for the 
corrections to scaling. The position ajmax — —5.4(1) of the 
maximum is insensitive with respect to these choices for the 
form of the corrections. In spite of this caveat the compari- 
son with the experimental data is encouraging. Note that the 
actual phase transition of the film occurs at a nonzero value 
of the bulk field. 



creases. For (-I--I-) [(+—)] BC the minimum [maximum] 
occurs above [below] Tc for all d. The shapes of the scal- 
ing functions in d = 2 and d = 3 exhibit an interesting 
resemblance. 

As we pointed out above, in the case of (H — ) BC the 
fluctuations of the order parameter are enhanced by the 
presence of a strongly fluctuating interface in the middle 
of the film. This results in a critical Casimir force which 
is generally stronger than in the case of (++) BC, for 
which there is no such an interface. This is reflected by 
the fact that the amplitude i^^ is larger than that of 

e.g., z?^"f''V|i??+"^l - 3.8 for the data sets ob- 
tained by accounting for the corrections to scaling ac- 
cording to Eq. (Uni) (case (i)) and Eq. ^ (case (ii)). 
Even though field-theoretical MFT per se does not pro- 
vide quantitative predictions for the overall amplitudes 
of the scaling functions , it yields the relation [47| 



-4t? 



(MFT) 

+ + 



i~2x) 



(27) 



and therefore predicts -d^™,^^^ /\'d^^'^^ \ = 4 and that the 

maximum of "d^ [minimum of z?++] occurs below [above] 

Tc- Thus MFT captures already quite well the qualita- 
tive and quantitative differences due to the presence or 
absence of an interface in the film. In addition, the fluctu- 
ations of such an interface, occurring in particular at low 

temperatures, cause the scaling function to decay to 

zero for x — > —oo more slowly than the scaling function 
1)++, which is clearly visible by comparing Figs.[9]and[T0l 



2. Dirichlet- Dirichlet boundary conditions 

In Fig. [TT] we show the MC data corresponding to 
g{y; L,2L, A) (see Eq. (fT2|) ) for the Ising model with 
(O, O) BC, realized by free surface spins. The L- 
depcndence of these data is quite pronounced and resem- 
bles that for the XY model with the same BC (compare 
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Fig. Ula))- On the other hand, the aspect ratio depen- 
dence appears to be relevant only in the narrow interval 
^2 < y < — 1 (see Fig. [T2|) . which is similar to the case 
of the XY model with periodic BC (compare Fig. [3]). As 
anticipated in Subsec. IIV Al the Ising model in a 3D film 
with Dirichlet-Dirichlet or periodic BC displays its 2D 
critical behavior at a critical point which is located on 
the bulk coexistence line _ff = at a size-dependent tem- 
perature TciL) such that Tc(L ^ oo) = Tc{l + ycL~'^''') 
[2], where j/c is a non-universal constant which depends, 
inter alia, on the BC. From extrapolating the MC data 
for Tc(L) reported in Table II of Ref. [Hi to L ^ oo one 
infers Uc.oo = —2.5(5) for the Ising model with (O, O) 
BC. As in the case of the XY model, the residual de- 
pendence on p observed in Fig. [12] might be due to the 
influence of the 2D phase transition for y ~ VcOO- Such 
a dependence cannot be captured by ansatze such as 
the ones considered so far, which assume that the cor- 
rections to scaling due to p 7^ are independent of x. 
Therefore, in order to achieve a good collapse of the data 
sets corresponding to different lattice sizes we account for 
corrections to scaling by following the procedure applied 
to the XY model with (O, O) BC, but we do not con- 
sider an aspect ratio dependence, i.e., we use the ansatze 
in Eqs. ^ (case (i)), and ^ (case (ii)) with 

ri_2 = 0. As a result of the fitting procedure in the inter- 
val X G [—7, —4] we find gi = 6.55(8) and = 2.35(3) in 
case (i), and 32 = —2.877(15) and g^ = 2.35(3) in case 
(ii) . Figure [13] shows the corresponding resulting esti- 
mates of the scaling function d{x) of the critical Casimir 
force with an excellent data collapse. As before, we find 
that i?(x) is affected by the choice of the functional form 
of corrections to scaling. In the two cases (i) and (ii) one 
finds estimates of i?(a;) which have the same shape but 
the overall amplitude is reduced by a factor R ~ 0.866 in 
case (ii) compared with case (i). 

Due to the residual dependence on the aspect ratio 
P^ a;„iin(p) and i?rnin(p) decrease upon decreasing p and 
therefore the values of iJmin and Xmin quoted above over- 
estimate the actual 'd-axmi.P = 0) and '(?mm(p = 0). The 
accuracy of our data does not allow us to study in more 
detail the Casimir amplitude Ao,o = *^(0)/2 (as we did 

for A_|_+ and in Fig. [8]), which turns out to be very 

small for (O, O) BC. Indeed the estimate from the par- 
tially resummed e-expansion is Ao,o = —0.0164 23 1, 
whereas MC simulations yield Ao,o = -0.0114(20) [ia j. 
However, from our data for the scaling function we can 
estimate ^0,0 = —0.014(8). The corrections of form (i) 
yield for the pronounced minimum of the scaling func- 



tion = -5.74(2) and d 



(i) _ 
mill 



i9(: 



whereas those of form (ii) result in x 



(ii) 



= -1.629(3) 
-5.73(4) and 



1? 



(ii) 



d{x 



.(ii) 



-1.41(1). 



In 2D the scaling functions obey the relation "doo (a^) = 
!?++(- x) [l^l. We note that in 3D this relation holds ap- 
proximately for the positions of the minima of the scal- 
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Monte Carlo data for 



goo{y = t(L - 
L,2L,A = {L/pf) (see Eq. (|12}, r = {T~Tc)/Tc) in 
the three-dimensional Ising model with (O, O) BC for L = 8, 
12, 16, 20, and for a fixed aspect ratio p — 1/6. The 2D crit- 
ical point of the film is located at y — yc,oo = —2.5(5), as 
inferred from extrapolating the data in Table II of Ref. [i^ 
to L 00. 
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FIG. 12: Monte Carlo data for goo{y = r{L - 
i)!/"; L, 2L, A = (L/pf) (see Eq. ([HJ , r = (T - in 
the three-dimensional Ising model with (O, O) BC for L = 12 
and various values of the inverse aspect ratio 1/p = yTKjL. 



than the scaling function for the (++) BC. For com- 
parison in Fig. [13] we provide the exact result for the 
two-dimensional Ising model (dashed line). In the in- 
set we show our MC data corresponding to the case 
(i) together with the scaling function obtained by us- 
ing the e-expansion in Ref. We note that it yields 
z9(0)/2 = Ao,o - -0.0118, which is larger than the es- 
timate Ao.o ~ —0.015 given in the same paper Q, and 
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obtained from dimensional interpolation; the latter value 
is still larger than the more recent theoretical estimate 
Ao,o - -0.0164 in Ref. [13]. 

In the case of Dirichlet-Dirichlct boundary conditions 
discussed here, the film exhibits the 2D critical behav- 
ior at T = Tc{L), corresponding to a universal value 
X* = yc{£.o)~^^'^ of the scaling variable x. Close to the 
temperature Tc{L), the free energy of the film is expected 
to exhibit the singularity - \T -Tc{L)\'^~°'^° , where a2D 
is the critical exponent of the specific heat of the two- 
dimensional system. This implies Q that the scaling 
function 'doo i^) of the Casimir force displays a singular- 
ity ~ \x-x*\'^-"^° a.tx = x*, i.e., ~ (a; - a;*)^ In |a; - a;*| 
for the Ising model. This singularity is too weak to be de- 
tectable by the present MC data. In Fig. [13] the gray bar 
indicates the value of Xqq = yc,oo{^o)~^^'^ = —7.6(1.3) 
and the associated uncertainty. Accordingly, the singu- 
larity is expected to occur on the left side of the pro- 
nounced dip. 

So far there are no experimental data available that 
would correspond to the Ising universality class with 
(O, O) BC. For experiments with binary liquid mixtures 
the (O, O) BC would correspond to walls which have no 
adsorption preferences, i.e., both components of the mix- 
ture are attracted equally by each surface. Effectively, in 
the limit of large film thicknesses, this can be achieved 
by chemically decorating the confining walls by stripes 
of equal width and alternating preferences for the two 
species of the binary liquid mixture (see Fig. 6 in Ref. [i^ 
for S = 1, for which within MFT the effective Casimir 
amplitude vanishes, corresponding to vanishing surface 
fields within MFT so that (O, O) BC hold). 



3. Periodic boundary conditions 

In the case of periodic BC the aspect ratio dependence 
of the Monte Carlo data for the Ising model (as for the 
XY model discussed in Subsec. flV Ap turns out to be rel- 
evant only in the vicinity of the minimum of the function 
gp{y', L, 2L, A) which is associated with the finite-size ef- 
fects close to the actual critical point of the thin film. 
Extrapolating the data in Table I of Ref. [1^ to L ^ cx) 
one infers that the shifted critical point corresponds to 
y = yc,p = —1.60(2). The fact that this type of finite- 
size dependence does not occur for the Ising model with 
fixed BC (see Fig. [7]) might be related to the different 
phase behavior below Tc in the latter case. For (++) BC 
the critical point is shifted off the bulk coexistence line 
_ff = to some value {Tc{L), Hc{L)) [HO] and hence in the 
vicinity of the minimum of the function (y; L, 2L, A) 
the corresponding bulk correlation length is smaller than 
the characteristic transverse length L II = ^/A. As already 
mentioned earlier, for (-1 — ) BC below Tc (but above the 
temperature of unbinding of this interface from one or 
the other surface) there exists a single film phase char- 
acterized by the OP profile displaying an interfacelike 
structure centered at the middle of the film [i^, [5l| . In 
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FIG. 13: Scaling function i9oo of the Casimir force for the 
three-dimensional Ising model with (O, O) BC and zero bulk 
field. The MC data refer to lattices with L = 8, 12, 16, 20 
and with a fixed inverse aspect ratio 1/p = 6. Corrections 
to scaling have been accounted for according to two different 
ansatze, provided by Eq. ((20} and Eq. (|2ip . and the corre- 
sponding numerical results are denoted by (i) and (ii), re- 
spectively. With corrections of the form (ii) , the shape of the 
resulting scaling function is almost indistinguishable from the 
one obtained with corrections of the form (i), but its overall 
amplitude is reduced by a factor R ~ 0.866. For comparison 
we show the exact result for the 2D Ising model (dashed 
line) and the mean-field prediction [2^ (dash-dotted line) nor- 
malized such that it yields the same depth of the minimum as 
the one of the MC data (i). In the inset we compare the MC 
data corresponding to the case (i) with the scaling function 
obtained from the e-expansion 01 . The gray bar indicates the 
value Xqq = —7.6(1.3) (and its uncertainty) of the scaling 
variable x corresponding to the occurrence of the shifted crit- 
ical point, inferred from extrapolating the data in Table II of 
Ref. in to L -> oo. 



this film phase the parallel correlation function f|| gov- 
erning the exponential decay of correlations along the in- 
terface is very large even for temperatures further away 
from Tc, i.e., ^|| ~ exp(L||/(4^)) with L|| = = Ly. i\\ 
gives rise to the aspect ratio dependence of the function 
5p(y;L,2L,A) for y < -4. 

In order to account for the corrections to scaling we fol- 
low the same procedure which we used for the XY model 
with periodic BC, i.e., we assume their L-dependence to 
be captured by Eq. (|22p at least within the range of sizes 
we are interested in. Accordingly, we focus on the data 
for the critical Casimir amplitude Ap and we fit them 
according to Eq. ([Ml) (case (iii), Eq. ^F^). The best 
fit parameters, based on all data points, are given by 
53 = 16.10(55) and Woff = 2.664(27) and the resulting 
curve is provided as a solid line in Fig. [Mj The associ- 
ated estimate for the asymptotic value Ap(i — > cx)) = 
Ap = -0.1520(2) agrees very weh with the MC resuh 
-0.1526(10) from Refs. 
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FIG. 14: MC data for the critical Casimir amplitude Ap{L) 
for the three-dimensional Ising model with periodic BC, as a 
function of the inverse film thickness L (on lattices with fixed 
inverse aspect ratio 1/p = 6). Due to L-dependent corrections 
to scaling, Ap depends on L and reaches its asymptotic value 
in the limit L ^ oo. The solid line corresponds to the best fit 
obtained by using the fitting ansatz given in Eq. ([22} in the 
interval < 1/L < 0.25. Our estimate (■) for the asymptotic 
value of the Casimir amplitude Ap(oo) compares very well 
with the previous MC result (•) from Ref. [23j . 

The scaling function iDp can now be determined by as- 
suming that Eq. (|19p . with r2 = and the parameters 
and Weff obtained from the analysis of Ap(L), effectively 
describes its corrections to scaling (case (iii), Eq. (HU), 
which actually leads to a very good data collapse in a 
wide range of temperatures. It also turns out that no 
corrections to the scaling variable x (see Eq. ^TE\\ ) are 
required in order to achieve it, i.e., ri, ~ 0. (Note, 
however, that corrections due to p ^ might be par- 
ticularly relevant within a certain range of the scaling 
variable x, see below.) The resulting scaling function iJp 
is presented in Fig. [TS] and it is based on a larger set of 
geometries of the simulation cell and with a better accu- 
racy than in our earlier work [l^ . The scaling function is 
in very good agreement with its previous determination 
in Ref. [la| based on the computation of the lattice stress 
tensor. The slight discrepancies might be due to the un- 
certainty in the normalization factor which had to be 
used in Ref. [ll] (see also Subsec HVAp . This agreement 
provides additional support concerning the reliability of 
our approach. 

Figurc[T5]prcsents also the comparison with the analyt- 
ical prediction of the recently proposed field-theoretical 
(FT) expansion up to 0(e^/^) [2l| (dash-dotted line) for 
X > 0. This latter prediction is now in better agree- 
ment with the MC data than the previous 0(e) field- 
theoretical result in Ref. 01 but still misses the onset of 
the formation of the minimum. (Figure 5 in Ref. (T^ 
compares the MC data with the 0(e) results, revealing 
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FIG. 15; Scaling function &p{x) of the critical Casimir force in 
the three-dimensional Ising model with periodic BC and zero 
bulk field. The data points refer to lattices with the inverse 
aspect ratio 1/p = 6. The corrections to scaling are taken into 
account according to Eqs. (jlQp and H22|) (case (iii)), see the 
main text). For comparison we show also the data set corre- 
sponding to the lattice with thickness L = 20 as investigated 
in Ref. the analytical prediction of Ref. [2l[ (dash-dotted 
line) for a; > 0, and results for 2D Ising model (dashed line) 
that we liave obtained numerically by using the transfer ma- 
trix method. Due to the self-duality of the 2D Ising model 
one has ■dp{—x) = &p{x) for d = 2 which allows for the oc- 
currence of two symmetric minima [SJ . We note that MFT 
yields ^p{x) = (solid line). The gray vertical line indicates 
the universal value x*p — —1.60(2) of the scaling variable x 
corresponding to the occurrence of the shifted critical point, 
inferred from extrapolating the data in Table I of Ref. [431 
L —> oo. 



a significant discrepancy for < a; < 4.) The estimated 

value of /P\0) = -0.39 from Refs. 0, [HI does not 
agree with our MC estimate t9p(0) = -0.3040(4). For 
the minimum of the scaling function we find the esti- 
mates Xmin = -0.681(1), l^inin = ^p{Xnun) = -0.329(1). 

Note, however, that for x ~ x^m the corrections due 
to p ^ are expected to be relevant. In order to sub- 
stantiate this statement we have determined the func- 
tion gp{y;L,2L,A = {L/pY) (see Eq. also from a 
set of data for lattices of thickness L = 10 and 15 with 
an inverse aspect ratio p^^ = 14, which can be com- 
pared with the corresponding data set from Fig. 1151 for 
which p~^ = 6. This comparison is presented in Fig. \W\ 
and clearly shows that, while the function gp is actually 
only slightly dependent on the thickness L of the lattice, 
there is a dependence on p which, however, is relevant 
only very close to Xmin- As mentioned above for the 
case of (O, O) BC, this latter dependence on p cannot be 
captured by ansatze such as the ones considered so far 
because they assume x-independent corrections due to 
p ^ 0- Therefore, similar to the case of {O, O) BC, due 
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FIG. 16: Aspect-ratio dependence of the function gp{y — 
t{L~ ^)^''' -1,21, a = {L/pf) (see Eq. ^) for the three- 
dimensional Ising model with periodic BC. For a fixed value of 
p, this function depends only weakly on L. By changing p, the 
function gp is affected mainly in the region —0.6 J/ J; 0. The 
2D critical point of the film is located a.t y = yc,p = —0.52(2), 
as inferred from extrapolating the data in Table I of Ref. [4g| 
to L ^ oo. 
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FIG. 17: Plot of the function gp{y\L,2L,A = {L/pf) (see 
Eq. p2p ') and the associated scaling function 6{y) (i.e., the lat- 
tice estimate of 6{y) = ^p{y / which is calculated by 
solving Eq. (|12|l ') iteratively for the three-dimensional Ising 
model with periodic BC. The data points refer to a lattice 
with L = 10 and 1/p = 6. The data have not been corrected 
for the corrections to scaling. The pronounced shoulder orig- 
inally present in gp is smoothed out upon calculating the 
associated scaling function. 



to the residual dependence on p, Xmiu{p) and t?min(p) de- 
crease upon decreasing p and therefore the values of i?niin 
and Xniin quoted above overestimate the actual values of 

l9min(/0 = 0) and Xmin. 

As in the case of Dirichlct-Dirichlet boundary condi- 
tions, the point at which the film exhibits the 2D critical 
behavior is located on the bulk coexistence line and cor- 
responds to a value x* of the scaling variable x, at which 
the scaling function is expected to display the weak sin- 
gularity ~ [x — x*)"^ \i\\x — X *\ (see Sec. |IVB2| above). In 
Fig. [12] the gray vertical line indicates the corresponding 
universal value x*p = yc,p{£,Q)~^^'' = —1.60(2). Accord- 
ingly, also in this case the singularity is expected to occur 
on the left side of the pronounced dip but cannot be de- 
tected by the present MC data. 

As a final remark we point out that for the Ising model 
with periodic BC the function gp{y; L,2L, A) exhibits a 
somewhat peculiar shape near Tc with a characteristic 
"shoulder" formed above the critical temperature. The 
procedure for retrieving the scaling function 9{y) [the 
lattice estimate of e{y) = 'dpiy / {^oY^'')] via Eq. ^ 
involves rescaling of the argument of the scaling function 
and removes the "shoulder" structure from the curve. 
The formation of this "shoulder" is related to the partic- 
ular shape of 'dp which on the left side of the minimum 
increases more steeply than on the right side of it (see 
Fig. [HI). 



V. SUMMARY AND CONCLUSIONS 

A. Summary 

We have presented important details of a novel general 
approach |12| to determine the universal scaling functions 
?9 of critical Casimir forces via MC simulations. We have 
applied this method (see Subsects. IIII Al and llll Bi as well 
as Fig. [T]) in order to study the scaling functions cor- 
responding to the three-dimensional Ising and XY bulk 
universality classes for a variety of universal boundary 
conditions in film geometries with varying thickness L. 
Corrections to scaling appear to be quite relevant in the 
range of sizes L we have investigated, which are strongly 
limited by the steeply increasing computational costs re- 
quired for larger systems. In spite of these difficulties, it 
is possible to analyze the corresponding MC data by as- 
suming suitable ansatze for corrections to scaling. Even 
if the final numerical determinations of the scaling func- 
tions are biased by these assumptions, they turn out to 
be consistent with the results of different numerical and 
analytical approaches and with all available experimental 
data. 

Our main results are the following: 

(1) We have obtained the Casimir scaling function 
i^oo for the three-dimensional XY model with (O, O) 
BC [(Dirichlet, Dirichlet) BC] (Fig. [J). Corrections to 
scaling have been accounted for by using two different 
ansatze, provided by Eq. (case (i)) and Eq. (PT]) 

(case (ii)). These choices of the functional form of cor- 
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rcctions to scaling have been dictated by the pronounced 
dependences on L and on the aspect ratio p of the sim- 
ulation cell which occur for this type of BC (Fig. El^a))- 
Both ansatze lead to a very good data collapse but the 
overall amplitude of the scaling function is reduced by a 
factor R ~ 0.9 in case (ii) compared to case (i). Our MC 
data compare very well with the corresponding experi- 
mental data for ^He films from Rcf. [1] and with the MC 
data of Ref. • For comparison also mean field results 
are provided. 

(2) The Casimir scaling function dp and the critical 
Casimir amplitude Ap have been obtained for the three- 
dimensional XY model with periodic BC (Figs. [6] and [5]). 
In this case, judged by the behavior of the generating 
function g introduced in Eq. (fT2|l . corrections to scaling 
are much less pronounced than in the case of {0,0) BC 
(Fig. [DJb)) and the aspect ratio dependence is relevant 
only in the restricted range of the scaling variable near 
the minimum of the scaling function (Fig. [3]). A very 
good data collapse is achieved by using the ansatz with 
the effective exponent o^cff = 2.59(4) (Eq. (case (iii)) 
and by neglecting the corrections to scaling due to the 
aspect ratio dependence (ri,2 = in Eqs. pS]) and p9|)). 
The shape of our MC data agree very well with the corre- 
sponding MC data of Ref. [1^ which, however, have left 
the amplitude undetermined. Our estimate for the criti- 
cal Casimir amplitude is Ap = —0.2993(7). By extend- 
ing the line of arguments of Ref. [11] to the present case, 
we have theoretically predicted the value •d^p^\—oo) = 
— C(3)/7r ~ —0.38 [see Eq. ([25)) ] at which the scahng func- 
tion i9p(x) saturates for a; ^ — oo. This value is con- 
firmed by the corresponding estimate —0.383(4) based 
on our MC data. 



(3) We have obtained the scaling functions -d^ 



^4 



and the corresponding Casimir amplitudes A-|__|-, A-j of 

the critical Casimir force in the three-dimensional Ising 
model with (++) and (H — ) BC, respectively, applicable 
for classical fluids (Figs. [HI [10] and El) . We find that in 
the critical regime the numerical data are practically in- 
dependent of the aspect ratio p (Fig.[7| but L-dependent 
corrections to scaling are rather important (Fig. [5]) . The 
presented scaling functions and Casimir amplitudes have 
been obtained by accounting for corrections to scaling 
according to Eq. ([^ (case (i)), Eq. ([^ (case (ii)), and 
Eq. ([^ (case (iv)) with ri,2 = (thus neglecting the 
dependence of the data on p). The final estimate of the 
scaling function is biased by the functional form assumed 
for the corrections to scaling; all considered cases provide 
a very good data collapse. The fitting ansatz in Eq. ([^^ 
describes very well the data for the Casimir amplitudes 
A^+/^ as a function of the film thickness L. Our es- 
timates for the asymptotic values of the Casimir ampli- 
tudes are A++ = -0.376(29) and A+_ = 2.71(2) which 
compare reasonably well with previous MC results from 
Ref. [2^ and with the results from the de Gennes-Fisher 
local- functional approach [l^ . Our results for the case of 
(H — ) BC compare well with recent X-ray scattering data 
for critical films of a classical binary liquid mixture Q. 



Moreover the MC data for the scaling functions -(9++ and 

have been used to calculate, within the Derjaguin 

approximation, the corresponding scaling functions for 
the critical Casimir potentials for the experimentally rele- 
vant geometry of a sphere near a planar substrate. These 
numerical results agree remarkably well with the experi- 
mental data for colloidal particles immersed in a critical 
solvent and close to a container wall [Tl| . 

(4) We have obtained the Casimir scahng function z?oo 
for the three-dimensional Ising model with (0,0) BC 
(Fig. [n|). For these BC the L-dependence of the MC 
simulation data is quite pronounced fFig. [TT|). similarly 
to the case of the XY model with the same (O, O) BC. 
The dependence on the aspect ratio is relevant only in the 
small range of the scaling variable near the minimum of 
the scaling function fFig. fT^ . Our data do not allow us to 
obtain a quantitatively accurate estimate of the Casimir 
amplitude, because of its very small value. Corrections to 
scaling have been accounted for according to the ansatze 
provided by Eq. ([20]) (case (i)) and Eq. (case (n)) 
with ri 2 = (thus neglecting the dependence of the data 
on the aspect ratio p). 

(5) The scaling function dp{x) and the critical 
Casimir amplitude Ap have been obtained for the three- 
dimensional Ising model with periodic BC (Figs. [TSl and 
rnj) . As in the case of the XY model with periodic BC, 
the aspect ratio dependence of the MC data appears to 
be pronounced only near the actual critical point of the 
thin film fFig. [TBI) . Therefore, the corrections to scaling 
have been accounted for in the same way as for the XY 
model with periodic BC, i.e., according to Eqs. and 
([22)) (case (iii)). The best fit for the L-dependence of the 
Casimir amplitude Ap has been obtained by using the 
ansatz given in Eq. (|26p . Our improved estimate for the 
value of the Casimir amplitude Ap = —0.1520(2 ) ag rees 
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very well with the previous MC result from Refs. 
The particular shape of the scaling function 'dp{x) around 
its minimum is rcfiected in the formation of a character- 
istic "shoulder" in the corresponding generating function 
gp above the critical temperature (Fig. [TT)) . 

B. Conclusions and outlook 

Our approach can be applied in order to study other 
experimentally relevant geometrical settings as well as 
the effect of chemically or geometrically inhomogeneous 
confining surfaces on the critical Casimir force. In the 
latter cases, even lateral critical Casimir forces are ex- 
pected to act in addition to the normal Casimir force in- 
vestigated here. This lateral force has been theoretically 
investigated for chemically [s^l and topographically 55 1 
patterned surfaces, whereas it has been experimentally 
studied for colloidal particles exposed to chemically pat- 
terned surfaces [56| . 

In addition to appliying our quantitative method to 
these cases, it is also desirable to perform more extensive 
and larger scale MC simulations in order to identify the 
origins of the corrections to scaling and to characterize 
them more accurately, possibly to the extent which is by 
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now achieved for bulk critical phenomena. This valuable 
knowledge would therefore allow an unbiased and thus 
even more accurate determination of the scaling functions 
of the critical Casimir force beyond the results presented 
here. Finally, beyond the application to thin *He films 
near the superfluid-normal fluid transition, our results for 
the three-dimensional XY model with (O, O) BC could 
be relevant for critical Casimir forces acting on Bose- 
Einstein condensates 1571. 
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APPENDIX A: CORRECTIONS TO SCALING 
AND FITTING PROCEDURES. 

In this appendix we describe the general strategy we 
have used in order to obtain the best fitted values of the 
parameters which control the corrections to scaling. The 
main problems one faces are to quantify the quality of 
a certain data collapse and then to choose the parame- 
ters which influence it in such a way as to optimize this 
quality. The estimation of the parameters and of the 
associated confidence interval proceeds as in the case of 
least-square fits with chi-square tests of the quality of the 
fit, but with the additional complication that the fitting 
function itself is not known and has to be estimated from 
the numerical data itself. 

In what follows we describe the procedure we have 
used in order to determine the best fit parameters which 
control the L-dependent corrections to scaling. On the 
same footing we have also treated the corrections due 
to a nonzero aspect ratio p ^ (see Subsec. IIIIPp . In 
full generality, assume that one seeks to determine, e.g., 
via MC simulations, the finite-size scaling function /i of a 
quantity ip which, in the absence of corrections to scaling, 
is expected to be a function of a scaling variable x only 
(which involves a suitable combination of temperature 
and size L of the system) so that 

%l;{x,L ^ oo) ^ h{x). (Al) 

For the time being we omit possible algebraic L- 
dependent prefactors of h. In the MC simulations one 
considers a set of N lattices of sizes Li, L2, 
and by varying the temperature one collects for each 
size Lk a discrete set of numerical values of ip with 



i = 1, 2 . . . which correspond to values Xk.j of the 

scaling variable x in the interval [xmin, a^max]- In this 
process the statistical uncertainty Aipkj associated with 
ijjk.j is also determined. From these quantities i/j^.j one 
intends to determine h, taking into account the presence 
of corrections to scaling. Due to them, ip is actually not 
a function of x only, but also of the size L of the sys- 
tem. In order to cope with this one therefore assumes 
the following functional structure: 

t^ix; L) - /i(L; h)hiML; h)x) (A2) 

where fi{L;ti) and f2{L\t2) capture the effects of 
the correction to scaling on the quantity "0 itself and 
on the scaling variable x, respectively. These func- 
tions depend on the size L of the system and on 
certain parameters ti , t2 which one would like to de- 
termine in such a way as to achieve the best data 
collapse for the function /i, obtained from the set 
of data points (/2(ifc; ^2)3::^^ , [./i(^fc, ii)]"Vfcj) =: 
iyk,jiti,t2),hkj(ti,t2)) for the various values of j and 
k, and as to take also into account the statistical er- 
ror Ahk.j(ti) := [fi{Lk,ti)]~^ Ai/jk,j associated with 

For each value Lfc we have interpolated the data set 
{xkjTipkj) in the interval [xmin , a^max] by using a cubic 
spline approximation. This way we have constructed a 
function ipk{x) with x S [smin, a^max] and with ipkixkj) = 
tpk.j- From this function we have calculated the corre- 
sponding Lfc-depcndent estimate hk{x;ti,t2) of h, given 
by 

hkivMM) = /i(ife;ii)"Vfc(y/2(ife;i2)"'), (A3) 

which fulfills (ii, i2); ti, ^2) = hk^j{ti,t2). In order 

to assess the quality of the data collapse and the quality 
of the fit we have actually to specify the function with 
which we would like to fit the data, which is the yet 
unknown scaling function h. In order to achieve this, we 
define an expected model function /lexpect as the average 
of the various hk, 

1 ^ 

^expect(y;ti,t2) = ^ /ifc (?/; ^1 , ^2) , (A4) 

fe=l 

which will then be fitted to the observed MC values by 
adjusting the parameters ti and t2- 

Accordingly, we calculate the "x^(^ii^2)" associated 
with the fitting of the data points (y^j (^i, ^2), /^fej (^1,^2)) 
with the function ft.cxpoct(j/; ^i, ^2): 

[^fe,j(ii,i2) — h 

expect {VkAtiMMM)? (A5) 
hh [^hkAt.)? 

Due to the non-trivial and non-linear dependence of the 
fitted data (and of the fitting function) on the parame- 
ters ti we cannot assume this quantity to play the same 
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role as a in more standard fitting procedures in which 
only the fitting function depends on the parameters one 
wants to estimate. Nevertheless, wc have heuristically 
made this assumption, i.e., that x^(ti, ^2) plays the same 
role as a x^, in order to determine the best fit param- 
eters and the associated confidence intervals. Accord- 
ingly we have proceeded as usual by determining the 
optimal fit parameters ti and ^2 which minimize the 



value of x^: X^(*i,i2) = min{tj_t3} x^(t]_, t2)-_ In order 
to estimate the statistical uncertainty At^ of ti we have 
determined that region of the plane (ti,t2) for which 
X^ih^h) < X^{hM) + 2.3 [111. The projection of the 
resulting region (typically of the form of an ellipse) onto 
the axis ti gives 2Ati, so that the estimate for the pa- 
rameters is of the form ii ± Aii. 
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